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SUMMARY 


Solid-state-batteries (SSBs) present a promising technology for next-generation batteries 
due to their superior properties including increased energy density and safer electrolyte 
design. Traditional SSB architecture features a Li-metal anode, a solid-state composite 
cathode and a stiff ceramic electrolyte for conduction of ionic species. While an attractive 
alternative, commercialization of SSBs faces a series of chemo-mechanical issues across 
1ts constituents, namely growth-induced fracture of the solid-state electrolyte (SSE) due to 
Li-metal deposition, interphase formation and associated large volumetric deformations at 
the anode-SSE interface, in addition to mechanical degradation across the various phases in 
composite electrodes. Computational modeling of SSBs is still at its infancy and constitutes 
the central theme of this thesis. We develop here a series of numerical frameworks seeking 
to elucidate the complex electro-chemo-mechanical processes coupled with mechanical 
fracture governing the behavior of both composite cathodes and solid-state electrolytes 
with Li-metal deposition. 

Both cathode and anode architectures may consists of a composite of active particles 
surrounded by a ceramic SSE matrix. During cycling, active particles undergo large vol- 
umetric changes against the stiff SSE, which can lead to formation of high mechanical 
stresses. Developing models for these architectures requires then a coupled electro-chemo- 
mechanical understanding of three critical constituents, namely active particles, solid-state 
electrolyte and the shared interface between the two. Key to such understanding is the role 
of mechanical damage on electrode performance. Towards modeling the role of mechani- 
cal damage on performance, we propose a novel coupled chemo-mechanical interface ele- 


ment, analogous to cohesive elements used in fracture mechanics. The framework enables 
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for modeling galvanostatic charging of composite electrodes and captures the continuous 
evolution of mechanical stresses, interfacial damage, and non-uniform current distribution 
across particles in a microstructure. We specialize on a Li,CoO2 (LCO) - Li¡oGeP2S¡2 
(LGPS) composite cathode and study the evolution of interfacial damage under varying 
material and microstructural properties. Specifically, we model electrolyte compositions 
with varying SSE stiffness and active particle volumetric change to assess their effect on in- 
terfacial stresses, mechanical damage, and overall electrochemical response. Subsequently, 
we discuss how variations in microstructural composition alter the state of interfacial dam- 
age. Both packing effects and particle size distribution are studied to understand how these 
factors impact integrity of the interface and performance. 

Towards modeling the phenomenon of Li-metal filament growth, we formulate a con- 
tinuum electro-chemo-mechanical gradient theory which couples phase-field damage and 
electrochemical reactions within the solid-state electrolyte. The proposed framework is 
fully-coupled with electrodeposition impacting mechanical deformation, stress generation 
and subsequent SSE fracture. Conversely, electrodeposition kinetics are affected by me- 
chanical stresses through a thermodynamically-consistent driving force that distinguishes 
chemical, electrical, and mechanical contributions. Importantly, the theory captures the in- 
terplay between crack propagation and electrodeposition by tracking the extent of damage 
and electrodeposition across the SSE via distinct phase-field variables, such that filament 
growth is preceded by and confined to fractured regions within the SSE. An attractive fea- 
ture of such approach is its ability to simulate the nucleation, propagation and branching of 
cracks and Li-filaments in arbitrary orientations. 

While the framework is general in nature, a specialization towards modeling the growth 
of Li-metal filaments in an inorganic Li7zLa3Zr2O;2 (LLZO) electrolyte is demonstrated. 
We showcase the capacity to capture intergranular and transgranular crack and Li-filament 
growth mechanisms, both of which have been experimentally observed. In doing so, we 


elucidate the manner in which mechanics and fracture of the SSE impact electrodeposition 
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kinetics and Li-filament morphology. From a manufacturing standpoint, we additionally 
elucidate the role of mechanical boundary conditions (1.e. mechanical confinement) on the 
rate of crack propagation across the electrolyte versus electrodeposition of Li-metal within 
cracks. Under specific mechanical boundary conditions, we demonstrate the capacity of 
the framework to capture the experimentally observed phenomenon of crack fronts prop- 
agating ahead of Li-metal filaments, as cracks traverse the entire electrolyte in advance of 
Li-deposits. 

In conjunction with experimental efforts, the theoretical frameworks proposed in this 
thesis seek to elucidate the inherent electro-chemo-mechanical processes which govern me- 
chanical integrity and ultimately electrochemical performance of all-solid-state batteries in 
an effort to expedite their successful commercialization for future energy storage applica- 


tions. 
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CHAPTER 1 
INTRODUCTION 


1.1 Background 


Solid-state-batteries (SSBs) present a promising technology for next-generation batteries. 
In recent years, commercialization of SSBs has attracted significant research attention ow- 
ing to their promise of superior performance compared to conventional Lithium-ion bat- 
teries (LIBs). Transition to SSBs has significant advantages including enhanced current 
density, faster charging time, wider electrochemical window, and improved safety. Con- 
trary to conventional LIBs which make use of graphite intercalation anodes, incorporation 
of Li-metal anode designs enables for a ten-fold increase in volumetric capacity, leading 
to a significant improvement in energy density (3860 mAh/g) [1, 2, 3]. Owing to their 
wider electrochemical window (0-5V), SSBs also enable for coupling of Li-metal anode 
with high-voltage cathode materials as demonstrated by several recent studies [4, 5, 6]. 
From a safety perspective, replacement of toxic and flammable organic electrolytes with 
highly conductive (> 107? S/cm), mechanically superior inorganic solid-state electrolytes 
(SSEs) can alleviate concerns associated with failure of batteries due to short-circuit or ig- 
nition [7, 8]. A particularly attractive SSB architecture with high energy density features a 
lithium metal anode, a composite solid-state cathode and an inorganic SSE separating the 
electrodes [9]. Fig. 1.1 highlights these three particular areas of interest in an all-solid-state 
architecture as well as the mechanical modeling issues which must be addressed in order 
to understand and commercialize these systems. 

Traditionally, SSEs are categorized into two groups, pertaining to organic solid poly- 
mers and inorganic solids (oxide or sulfide based). Among other properties, successful 


integration of SSEs in SSBs requires meeting particular chemical standards, including neg- 


ligible electronic conductivity, high ionic conductivity, good electrochemical compatibility 
as well as several mechanical standards including appropriate stiffness and fracture tough- 
ness to preserve mechanical integrity. While liquid electrolytes traditionally have a higher 
ionic conductivity and allow for excellent wetting of electrode surface, comparable ionic 
conductivity values have been achieved for several ceramic SSEs as reported in multiple 
studies [10, 11, 12, 13]. Resolution of sluggish ionic conductivity, while a major improve- 
ment in research efforts to integrate SSEs in next-generation battery designs, addresses one 


of the multiple challenges hindering commercialization of SSBs. 
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Figure 1.1: Architecture of an all-solid-state battery including a Li-metal anode, a solid- 
state electrolyte (SSE), and a composite cathode. Key areas of investigation: (a) behavior 
of conversion electrodes undergoing sharp interface reactions; (b) evolution of composite 
cathodes including damage; (c) Li-filament growth within cracks; and (d) formation and 
evolution of interphases at the SSE/metal interface. Composite figure reproduced with 
permission from [9, 14, 15, 16, 17, 18, 19, 20]. 


A critical challenge here remains the ability to control chemo-mechanical instabilities 
and contact loss at the electrode-SSE interface [21, 22]. While research in understanding 
of solid/liquid interfaces is mature, the community’s understanding of solid/solid electro- 
chemical interfaces is still at its infancy. A multi-physics experimentally-guided under- 
standing of the complex nature of solid-solid interfaces and their evolution is thus critical 


to resolve a series of issues presently hindering commercialization of SSBs and consti- 


tutes the central theme of this thesis. In particular, as part of this thesis, we propose a 
series of novel theoretical frameworks seeking to model and elucidate the complex coupled 
electro-chemo-mechanical phenomena and associated degradation mechanisms in two pri- 
mary constituents i) composite electrodes and ii) solid-state electrolyte and its fracture due 


to Li-filament growth. 


1.1.1 Modeling of the electro-chemo-mechanical behavior of composite electrodes 


Composite cathodes, as illustrated in Fig. 1.1, present a critical constituent where me- 
chanics plays a key role. By construction, composite cathodes rely on a densely packed 
architecture of active material, surrounded by a conductive SSE to support ionic trans- 
port, in addition to complementary additives for electronic conduction [23, 24, 25]. The 
active material undergoes volumetric expansion/contraction during insertion/extraction of 
ionic species, which, due to the confined nature of the all-solid-state architecture, can lead 
to generation of high stresses [26, 27]. Batteries employing liquid electrolytes tradition- 
ally isolate volumetric strains at the active particle level, given the incapacitation of liquid 
electrolytes to transfer stresses. As a result, in current battery designs, any morphological 
changes of the active particles do not translate to interfacial delamination at the electrode- 
electrolyte interface, which continues to support ion transfer. 

The confined nature of solid-state electrodes allows for transmission of stresses through 
the stiff solid framework, resulting in a complex coupled chemo-mechanical behavior, 
which is prone to damage through failure of the active particles [28, 29, 30], the SSE 
itself [31], or delamination of the interface between the two [17, 32]. Fig 1.2(b-d) illus- 
trates the different degradation mechanisms across the constituents of a composite cathode 
for all-solid state battery applications, limiting in turn their electro-chemical performance 
as also experimentally reported [16, 32, 33]. 

In addition to the mechanical behavior, there are concurrent electrochemical processes 


as current (i.e. amount of charge) locally varies between particles in a composite electrode. 


From a modeling perspective, there is a need to capture multi-particle interactions—both 
mechanical and chemical —for these densely packed architectures in order to understand 
how variations in material properties and microstructure composition affect mechanical 
degradation and ultimately electrochemical performance. At present, modeling of the con- 
current chemical interactions between active particles arising from diffusion of Li-ions 
through the electrolyte remains largely unexplored, with most frameworks resorting to 
study of single-particle electrode domains [29, 30, 34, 35]. An inherent obstacle in mod- 
eling of chemical interactions between particles is implementation of true galvanostatic 
conditions. During galvanostatic charging, total current over all particles is held constant, 
but individual current distribution at each particle surface is unconstrained and depends on 
the local stress-coupled non-linear reaction kinetics. 
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Figure 1.2: (a) Architecture of a solid-state composite electrode consisting of a conglom- 
erate of active particles, solid-ionic conductor and binders for electronic conduction. Me- 
chanical fracture across constituents of a composite cathode, namely fracture at: (b) active 
particle bulk; (c) active particle-SSE interface; and (d) solid-state ionic conductor. Com- 
posite figure reproduced with permission from [17, 28, 31, 32, 36, 37] 


To address this limitation, in chapter 2 of this thesis, we formulate and numerically im- 
plement a theoretical framework based on application of chemo-mechanical interface ele- 
ments — analogous in form to conventional cohesive elements used in fracture mechanics — 
to model the stress-coupled non-linear kinetics of ionic transport at the particle/electrolyte 
interface. We showcase the capacities of the proposed surface elements to successfully 
capture galvanostatic charging conditions and cross-particle interactions by first modeling 
the electro-chemical performance of innovative LIBs electrode designs, in which interac- 
tions are purely of chemical nature. Subsequently, we deploy this new element in chapter 3 
to model galvanostatic charging of composite electrodes for all-solid-state battery applica- 
tions and study their electrochemical performance with variation in both chemo-mechanical 
properties as well as microstructure composition (i.e. active material infill, particle size dis- 
tribution). From a design standpoint, the proposed framework is particularly utilitarian and 
can help elucidate manufacturing efforts in terms of optimal microstructure composition to 


minimize damage and enhance electrochemical performance. 


1.1.2 Modeling of interphase formation at the electrode-solid electrolyte interface. 


Furthermore, unwanted electrochemical reactions between the Li-metal anode and SSEs 
during cyclic charge/discharge of the battery can produce compounds, which lead to forma- 
tion of a non-homogeneous interphase layer of low ionic conductivity. This in turn can sig- 
nificantly limit ion transfer kinetics and increase interfacial resistance as reported in several 
representative SSEs [5, 38, 39, 40]. Simultaneously, formation of non-homogeneous in- 
terphase regions is accompanied by volumetric expansion, which in the context of densely 
packed all-solid-state architectures can induce high stresses, resulting in fracture and degra- 
dation of the SSE [18, 41]. Occurring in the confined all-solid-state architecture, this vol- 
ume expansion introduces internal stresses, leading to possible fracture and degradation of 
the solid electrolyte. 


Tippens et al. [18] have reported the formation of cracks in a Liy,4Alo,4Ge1,6(PO.)3- 


(LAGP) solid-state electrolyte using X-ray tomography, interestingly arguing that it is the 
presence of these cracks and the subsequent loss of contact that cause high impedance in the 
battery operation, rather than the transport properties of the interphase itself. Fig. 1.3(a) il- 
lustrates the formation of the thermodynamically-unstable interphase at the Li-metal-LAGP 
interface, while a time sequence of the evolution in fracture of the SSE incurred by volu- 


metric expansion of the interphase is shown in Fig. 1.3(c). 
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Figure 1.3: (a) Formation of a thermodynamically-unstable interphase at the Li-metal 
anode-LAGP solid electrolyte interface. (b) Finite element model of the swollen inter- 
phase, showing (1) contours of reaction coordinate (left) and (11) evolution of radial and 
circumferential stresses across the LAGP pallet, leading to formation of cracks (right). 
(c) Time sequence of the fracture network induced by sequential volumetric growth of 
the interphase. (d) Direct correlation between extent of crack propagation and increase in 
impedance. Reproduced with permission from [18]. 


From a modeling perspective, the same team also studied the stress distribution in the 
aforementioned solid state cell using a diffusion-deformation continuum framework by Di 
Leo et al. [42] as shown in Fig. 1.3(b). It should be noted, however, that the model used 
to simulate the interphase formation and evolution relied on an ad-hoc description of the 
interphase, where its position was determined apriori rather than evolved according to the 


chemo-mechanical physics of the system. Similar volumetric changes incurred by forma- 
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tion of thermodynamically-unstable interphases at the Li-metal anode in contact with the 
solid electrolyte have been corroborated by Zhang and co-workers [43]. As reported in 
their work, the chemically formed interphase strongly affects the mechanical integrity of 
the LAGP solid electrolyte, pulverizing the LAGP pellet with repeated microcrack forma- 
tion on the side of the interphase. 

While extensively explored experimentally, at present, capturing the phenomenon of in- 
terphase formation and its evolution remains one of the challenges for SSE realization from 
a modeling perspective. Other continuum formulations in chemo-mechanics of solid inter- 
phases, such as the one by Rejovitzky et al. [44], are more suitable for modeling the growth 
of solid interphase in conventional LIBs, as opposed to structural and phase changes asso- 
ciated with interphase formation at the anode-electrolyte interface in SSBs; treatment of 
growth models in LIBs is fundamentally different from a reaction model relevant to SSBs. 
There are, however, models in the mechanics community from which useful formulations 
for interphase growth can be drawn, as they describe physical phenomena similar in nature 
to interface reactions in SSBs. Particularly, the research conducted on modeling oxidation 
in thermal barrier coatings [45, 46, 47], and modeling of hydrogen-assisted cracking in 
metals [48, 49] can serve as potential starting points for building models specific to inter- 
face reactions in SSBs. Nevertheless, at present, a diffusion-reaction-deformation theory 
capable of accommodating chemically-induced phase transformations is currently lacking 
in the mechanics literature. 

A key constituent in development of such frameworks is also the capacity to model the 
formation and propagation of cracks across the solid conductor induced by large interphase 
volumetric expansions. Phase-field damage formulations have been extensively reported 
in the literature with applications across a multitude of problems including active material 
cracking [36, 50, 51], stress-corrosion cracking [52, 53], hydrogen embrittlement [48, 54], 
soft-matter [55, 56, 57] etc. These models are easily transferable to modeling of SSE 


fracture induced by volumetric expansion of the interphase and importantly do not require 


prior definition of the crack path. Instead, cracks are free to evolve in any arbitrary direction 
driven by a thermodynamically-consistent fracture driving force, provided a critical energy 


condition for damage initiation is met. 


1.1.3 Modeling of Li-metal filament growth across solid-state electrolytes 


Among the numerous degradation mechanisms listed in Fig. 1.1, Li-metal filament growth 
upon fracture of the solid-state electrolyte constitutes one of the most intricate challenges in 
SSBs, resorting to the interplay of both mechanical and chemical forces. Contrary to den- 
drite growth processes in liquid electrolyte systems, growth of metal filaments across SSEs 
in solid-state architectures depends on a fundamentally different mechanism, hypothesized 
to be centered on the morphology of the SSE microstructure. Filament protrusions can ini- 
tiate at perturbations of the interface or microstructural heterogeneities (i.e. pores, cracks, 
grain boundaries) and subsequently grow through the SSE microstructure by fracturing the 
latter and creating fresh sites for continuous, non-uniform Li-metal deposition towards the 
cathode end. Fig. 1.4(a) illustrates the potential formation of surface imperfections in the 
form of voids at the Li-anode-SSE interface with continuous cycling of the solid-state bat- 
tery (1.e plating/stripping). Formation of voids at the Li-anode leads to current constrictions 
at the interface, localizing deposition at these sites. This in turn facilitates nucleation and 
subsequent propagation of Li-protrusions across the solid electrolyte with fracture of the 
latter. Eventually, as Li-metal protrusions form a full pathway from the anode towards the 
cathode, the battery short-circuits [20, 58, 59, 60]. It is thus critical to understand from both 
an experimental and modeling perspective the interplay between non-uniform electrodepo- 
sition at the Li-metal interface and heterogeneity of the SSE microstructure, in particular 
the size and distribution of surface defects, rate of lithium deposition at the Li/SSE inter- 


face, local stress state and the viscoplastic deformation of Li-metal [22, 61, 62, 63, 64]. 


Loading Direction 


Figure 1.4: (a) Formation of voids at the Li-metal anode interface with continuous plat- 
ing/stripping. Voids formation at the interface induces current constrictions at these sites, 
facilitating nucleation and subsequent propagation of Li-metal filaments with fracture of 
the solid-electrolyte. Reproduced with permission from [22]. (b) Transgranular growth 
of Li-metal filaments across a single-crystal LLZO microstructure. Reproduced with per- 
mission from [65]. (c) Preferential deposition of Li-metal along the grain boundaries of a 
polycrystalline LLZO microstructure. Reproduced with permission from [58]. 


From an experimental standpoint, filamentary protrusions across both compliant [20, 
66, 67, 68] and stiff oxide-based electrolytes [58, 59, 60, 69, 70, 71] have been extensively 
reported above a critical current density. These works have nullified the findings of Mon- 
roe and Newman that a region of stability for suppression of metal filaments growth exists, 
provided the shear modulus of the SSE is at least twice larger than that of Li-metal [72]. Ad- 
vanced characterization techniques including tomographic imaging and spectroscopy have 
concurrently revealed propagation of Li-metal filaments across SSE microstructures either 
1) transgranularly, or through cracks typical for single-crystal electrolytes [20, 60, 69, 73] 
and ii) integranularly or through the grain boundaries, typical for polycrystalline SSEs [58, 
59, 66, 74, 75] as shown in Fig. 1.4 (b-c). Recent studies have additionally hypothesized 
the potential for isolated Li-metal deposits to form within the bulk SSE microstructure due 
to the favorable SSE electronic properties and presence of trapped electrons. Upon reac- 
tion with Li-ions, these trapped electrons may produce isolated Li-metal deposits, seeding 


subsequent growth of metal filaments across the SSE microstructure [76, 77]. 


While an understanding of the complex interplay of electro-chemo-mechanical phe- 
nomena governing the growth of Li-metal filaments has been extensively explored from an 
experimental standpoint, development of a microstructure-resolved description of the pro- 
cess Of Li-metal filament growth through the SSE microstructure is lacking and necessary, 
with current frameworks resorting primarily to modeling growth of filaments as pressurized 
cracks under a linear-elastic fracture mechanics framework [78, 79, 80, 81, 82, 83]. While 
clearly important, these models are limited to modeling the conditions under which an ex- 
isting filament might propagate. They can not capture the dynamic evolution of Li-metal 
filaments. 

Only recently, mature continuum frameworks to model the electro-chemo-mechanical 
processes that govern non-uniform electrodeposition at a Li-metal/SSE interface in the 
presence of imperfections have been proposed. In [84], Narayan and Anand propose a con- 
tinuum framework to model non-uniform growth due to uneven deposition at the Li-anode 
coupled with elastic-viscoplastic behavior of Li-metal. Given the complexity of modeling 
freshly-deposited Li-layers from a finite element standpoint, the authors propose an analo- 
gous mechanical problem of swelling and de-swelling of a thin “interphase layer” at the Li- 
metal/SSE interface to reproduce the electro-chemical process of plating/stripping. How- 
ever, the framework is limited in its capacity to directly model the ongoing electrochemical 
reaction kinetics driving the growth of Li-metal filaments, and can not subsequently predict 
the morphology and evolution of the resulting Li-filaments. Alternate continuum frame- 
works, modeling Li-metal protrusions across the SSE as climbing dislocations have been 
recently proposed by Shishvan et al. [85, 86], importantly showing the capacity to predict 
the experimentally reported critical current density and growth rate of Li-deposits across 
SSEs consistent with experiments. However, at present, modeling frameworks accounting 
for the role of microstructure heterogeneity on both initiation and propagation of Li-metal 
filaments across SSEs are limited. 


Pertaining to initiation, only recently Zhao and co-workers importantly proposed a 
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phase-field formulation, which enables for modeling the evolution of the metal electrode- 
SSE interface with the potential for void formation during continuous plating/stripping 
[87]. Critically, formation of voids at the interface localizes reaction kinetics, leading to 
formation of current hot spots and potential nucleation of Li-metal filaments at these sites. 
Consistent with experiments, the authors predict the formation and evolution in void mor- 
phology during electrochemical cycling (1.e. plating/stripping) with applied stack pressure, 
and elucidate the interplay between a series of phenomena including bulk diffusion, Li 
dissolution/deposition, creep, and the nucleation and annihilation of vacancies. Similar 
treatments, investigating the formation and evolution of voids at the interface are reported 
in [88, 89]. 

Nevertheless, at present, an electro-chemo-mechanical framework which captures: 1) 
the nucleation and evolution of Li-metal filaments, ii) the concurrent fracture of the solid 
host, and 111) the coupling between mechanical deformation, stress and Li-metal electrode- 
position kinetics is missing. Of particular importance in such a framework is also the ability 
to model the role that heterogeneities in the SSE microstructure (i.e. pores, grain bound- 
aries) play on the coupled fracture behavior of the SSE and concurrent Li-metal filament 
growth due to electrodeposition. 

To address current modeling limitations, in chapter 4 of this thesis, we propose a 
thermodynamically-consistent, electro-chemo-mechanical gradient theory coupled with dam- 
age, which models the phenomenon of Li-metal filament growth inside SSEs. Critically, 
the proposed framework elucidates the role that heterogeneities in the SSE microstructure 
(i.e. pores, cracks, grain boundaries) play on the coupled fracture of the SSE and concur- 
rent Li-metal filament growth due to electrodeposition. In particular, through the proposed 
framework, we demonstrate the capacity to capture integranular and transgranular growth 
mechanisms of Li-metal filaments, both of which have been experimentally observed [19]. 
Altogether, commercialization of next-generation all-solid-state battery architectures relies 


on development of high-fidelity models across the various components of the cell, capable 
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of capturing the complex interplay between electro-chemo-mechanical phenomena coupled 
with fracture. Towards this goal, the research work outlined in this thesis seeks to provide a 
series of theoretical frameworks — which in conjunction with experimental efforts — can en- 
hance our understanding of the coupled phenomena governing instabilities in all-solid-state 
architectures in an effort to accelerate their successful commercialization. 

The overview on modeling aspects of all-solid-state batteries was adapted with permis- 


sion from the review work published in: 


e Bistri, Donald, Arman Afshar, and Claudio V. Di Leo. "Modeling the chemo-mecha- 
nical behavior of all-solid-state batteries: a review.” Meccanica 56 (2021): 1523- 


1554. 
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CHAPTER 2 
A NON-LINEAR KINETICS INTERFACE ELEMENT FOR MODELING OF 
MULTI-PARTICLE INTERACTIONS IN LI-ION BATTERY ELECTRODES 


Modeling of chemo-mechanical interactions between active particles in battery electrodes 
remains a largely unexplored research avenue. Of particular interest from a modeling stand- 
point is the ability to model the distribution in current (i.e. charge), which may vary across 
active particles under galvanostatic (1.e. fixed current) charging conditions. This in turn 
depends on the local stress-coupled electrochemical potential, and may also be affected 
by mechanical degradation. In this chapter, we formulate and numerically implement a 
constitutive framework, which captures the complex chemo-mechanical multi-particle in- 
teractions in electrode microstructures, including the potential for mechanical degradation. 
A novel chemo-mechanical surface element is developed to capture the local non-linear 
reaction kinetics and concurrent potential for mechanical degradation. We specialize the 
proposed element to model the electrochemical behavior of both (1) a novel anode de- 
sign for application in LIBs and (ii) a next generation all-solid-state composite cathode, 
where mechanical interactions are particularly important. In modeling these electrodes, 
we demonstrate the manner in which the proposed simulation capability may be used to 
determine optimized electro-chemical and mechanical properties as well as the layout of 
the electrode microstructure, with a focus on minimizing mechanical degradation and im- 
proving electrochemical performance. The work presented in this chapter is adapted with 


permission from and published in: 


e Bistri, Donald, and Claudio V. Di Leo. ”Modeling of chemo-mechanical multi- 
particle interactions in composite electrodes for liquid and solid-state Li-ion batter- 


1es.” Journal of The Electrochemical Society 168.3 (2021): 030515. 
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2.1 Introduction 


Lithium-ion batteries (LIBs) have become the dominant energy storage source in a number 
of areas including consumer electronics and electric vehicles [1]. A number of research 
avenues exist for the continued improvement of LIBs including improved charging times, 
increased operational lifespan, expanded capacity, and improved safety [90, 91]. In pursuit 
of these improvements, research efforts have shifted from conventional liquid LIBs towards 
next generation all-solid-state battery (SSB) architectures with enhanced current density (~ 
3860 mAh/g), wider electrochemical window (0-5 V), and improved safety [2, 3, 5, 6, 9]. 
However, battery designs, be it for conventional liquid LIBs or next-generation SSBs, suffer 
from inherent mechanical degradation across different cell constituents, which in turn limit 
their electrochemical performance. Upon intercalation/de-intercalation of lithium ions, cur- 
rent electrode designs undergo large inhomogeneous volumetric changes, which give rise 
to large mechanical stresses and inelastic effects [26, 27, 92, 93]. This in turn can lead to 
mechanical degradation, which for example may take the form of fracture of the individual 
constituents or debonding of active particles from their surrounding matrix, as has been re- 
ported in the literature [94, 95, 96, 97, 98]. The role of mechanical stresses and mechanical 
degradation on electrochemical performance is even more critical in all-solid-state archi- 
tectures owing to the presence of a stiff solid electrolyte. Aside from conventional fracture 
mechanisms at the electrode level [31, 99, 100], inherent mechanical instabilities associated 
with: i) metal filament growth [58, 60, 101], and ii) interphase formation and interphase 
induced fracture of the solid-state electrolyte [18, 41] have hindered commercialization of 
next generation SSBs as outlined in chapter 1. 

Research to date on the coupling between mechanical degradation and electrochemical 
performance remains mostly limited to study of single-particle systems (cf. [95, 102]). As 
such, our understanding of the role of mechanical degradation on electrochemical perfor- 


mance, including the effect of mechanics on interfacial kinetics, remains rather limited. 
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From both an experimental and modeling perspective, it is thus critical to understand the 
complex interplay of mechanical deformation and chemical processes. This is particularly 
critical in the context of densely-packed microstructures traditionally employed in design 
of composite electrodes for commercial batteries. 

Conventional composite electrode architectures build upon a compact conglomerate 
of active particles confined on their exterior by a matrix supporting percolation pathways 
for ionic conduction in addition to complementary additives for electronic conduction [16, 
103, 104]. Owing to their compact nature, composite electrodes experience a complex in- 
terplay of chemo-mechanical interactions as Li-ions are inserted/extracted from the active 
particles, thus leading to emergence of a complex stress field as the particles swell/deswell 
against the confining matrix. In addition to the mechanical behavior, there are concur- 
rent electrochemical processes as current locally varies between particles in a composite 
electrode. From a modeling perspective, there is a need to capture multi-particle interac- 
tions — both mechanical and chemical — for these densely packed architectures in order 
to understand how variations in material properties and microstructure composition affect 
mechanical degradation and ultimately electrochemical performance. 

Significant research efforts have been undertaken from a theoretical standpoint to un- 
derstand the chemo-mechanical behavior of single active particles, cf. Zhao et al. [105]. 
Modeling frameworks have been developed to different levels of complexity for both in- 
tercalation electrodes with deformations usually confined within the elastic regime [34, 35, 
106, 107, 108, 109], and conversion electrodes undergoing volumetric strains in the or- 
der of more than 100% and experiencing large elastic-plastic deformations [14, 42, 110]. 
Theories coupling large elastic-plastic deformations with large volumetric swelling due to 
diffusion of Li-ions have been formulated [93, 111, 112, 113] and specialized to capture 
the chemo-mechanical behavior of single Si-particles [114]. A key constituent is the abil- 
ity of these frameworks to capture the experimentally measured electrochemical response. 


Following a combination of experimental measurements and numerical simulations in a 
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Li-ion half cell, Bucci et al. [115] adopted the theory of Bower et al. [111] to characterize 
the behavior of amorphous thin film Si-electrodes. The effect of elastic-plastic deforma- 
tion on electrochemical performance of geometrically complex a-Si anodes experimentally 
realized by Wu et al. [116] was investigated in the work of Di Leo et al. [42]. 

While much has been done from a theoretical standpoint to model the chemo-mechanical 
behavior of single-particle systems, few studies have considered the combined multi-particle 
behavior of composite electrodes. In modeling the multi-particle behavior of electrodes, 
one must consider both mechanical and chemical interactions which occur between parti- 
cles. Mechanical interactions arise from transfer of forces between particles either through 
direct particle-to-particle contact or through contact with the surrounding matrix. Early 
modeling efforts aimed at capturing mechanical interactions between active particles in 
electrodes centered on applications to traditional liquid LIBs. Experimentally, Lee et 
al. [117] studied mechanical interactions in Si-nanopillars using both in-situ and ex-situ 
characterization techniques to explore how mechanical interactions among neighboring 
nanopillars impacted reaction kinetics and their susceptibility to fracture upon lithiation. 
Numerical finite-element simulations of a 2D porous graphite/binder electrode microstruc- 
ture were performed by Rahani et al. [118] to model mechanical interactions between 
active particles and PVDF binders. In a similar effort, Higa and Srinivasan [119] investi- 
gated the role of both particle size and binder stiffness on the mechanical stresses induced 
in active particles, reporting a decrease in stress with decrease in particle size and binder 
stiffness. Garcia et al. [120] modeled 2D porous electrodes of Li,Mn2O, spherical parti- 
cles under a linear-elastic finite-element framework. From a design standpoint, the effect 
of both packing density and particle size were studied. Microstructure configurations with 
dense particle aggregates were assessed to be detrimental to electrochemical performance, 
lowering the utilization of active material and leading to higher mechanical stresses induced 
in the particles. 


Recently, numerical simulations of 3D reconstructed composite electrode microstruc- 
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tures under a continuum framework have been reported by Zhao and co-workers [103, 104]. 
The interplay between mechanical stresses and Li-ion diffusion was explored to quantify 
its role on capacity for two representative microstructures including: 1) an NMC cathode 
undergoing small volumetric changes upon Li insertion, and 11) an SnO anode where volu- 
metric changes upon lithiation are more significant. These simulations were later expanded 
by the same group to model damage both at the particle level and particle interface in NMC 
composite cathodes [121, 122]. Efforts of a similar nature are also reported by Hoffman 
and co-workers [123], where stochastic microstructure models were used for representative 
graphite anodes. 

While these works represent a step towards modeling of complex chemo-mechanical 
interactions in composite electrodes, development of a general computational framework, 
along with a suitable numerical implementation, capable of capturing appropriate galvano- 
static charging conditions remains at its infancy and is addressed in this chapter. Of interest 
here is thus the development of a general computational framework, which captures the role 
of local variations in chemo-mechanical properties and microstructure composition on the 
mechanical integrity and electrochemical performance of different electrode microstruc- 
tures. 

Within the context of all-solid-state batteries, the most notable efforts in modeling of 
composite electrodes — with a focus on the role of mechanics — are by Bucci and co- 
workers [31, 99, 100]. In Bucci et al. [99], a joint analytical and finite-element framework is 
reported to model the effect of stresses in Silicon active particles within an SSE composite 
electrode. A critical Silicone to SSE volume fraction was reported above which mechan- 
ical stresses due to multi-particle interactions have a dominant effect on electrochemical 
performance and one is under-utilizing the added active material from a capacity stand- 
point. This work was subsequently extended to account for mechanical damage through 
use of cohesive elements both at the active particle-electrolyte interface [31] and the SSE 


itself [100]. Interestingly, compliant electrolytes were assessed to be more prone to micro- 
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cracking. Additionally, fracture initiation at the active particle-SSE interface was reported 
for most combinations of SSE material properties, provided the active particles undergo 
7.5% change in volume upon lithiation. From a design standpoint, the authors assessed 
initiation of fracture to occur regardless of the size of active particle, claiming that nanos- 
tructured electrodes need not perform better from a mechanical standpoint. In the context 
of all-solid-state batteries, the role of mechanical interactions on transport properties has 
been recently explored in the work of Behrou and Maute [124]. Here, the effect of damage 
on mechanical behavior and effective transport properties was studied using a non-local 
damage model with both mechanical stiffness and diffusivity of active particles penalized 
with damage evolution. In the same spirit, the role of mechanics on effective transport 
properties and electrochemical performance of SSBs with increase in active material was 
also investigated in the work of Al Siraj et al. [125]. A critical active material volume 
fraction was determined above which electrochemical performance ceases to improve as 
conduction pathways are significantly decreased. 

While research efforts have been devoted to capturing mechanical interactions at the 
electrode microstructure, modeling of the concurrent chemical interactions between active 
particles arising from diffusion of Li-ions through the electrolyte remains largely unex- 
plored. An inherent obstacle in modeling of chemical interactions between particles is im- 
plementation of true galvanostatic conditions. During galvanostatic charging, total current 
over all particles is held constant, but individual current distribution at each particle surface 
remains unconstrained and depends on the local stress-coupled non-linear reaction kinetics. 
This chapter is dedicated to formulating and numerically implementing a novel theoretical 
framework based on use of a chemo-mechanical interface element — analogous in form to 
conventional cohesive elements used in fracture mechanics — to model the stress-coupled 
non-linear kinetics of ionic transport at the particle/electrolyte interface. We deploy this 
new element to model galvanostatic charging of composite electrodes and study their elec- 


trochemical performance. The proposed interface element is introduced at the surface of 
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active particles within the numerically discretized composite electrode microstructures to 
capture the non-linear stress-coupled reaction kinetics. Further, mechanical degradation of 
the interface is concurrently modeled — accounting for a loss in reactivity due to dam- 
age — to assess the interplay of mechanical stresses and damage on the behavior of the 
composite electrode. 

We present the theoretical and numerical framework in Sect. 2.2. This section details 
the continuum description of a representative volume element (RVE) of the microstructure, 
Sect. 2.2.1, followed by a detailed description of the constitutive modeling and numerical 
discretization of particles interface in Sect. 2.2.2 and particle bulk in Sect. 2.2.3. Subse- 
quently, we showcase the capabilities of our numerical framework in Sect. 2.3 through a 
simple demonstration of a two-particle RVE. We then apply our numerical framework to 
model experimentally relevant electrode microstructures. First, in Sect. 2.4, we investigate 
a traditional liquid LIB electrode design composed of multiple double-walled hollow nan- 
otubes. Here, only chemical interactions between the multiple nanotubes are considered 
and we explore the role of material properties and interfacial damage on electrode perfor- 
mance. As later detailed, in chapter 3, we specialize the proposed numerical framework 
and model a representative solid-state composite cathode where mechanical load transfer 
through the stiff solid-electrolyte is critical and can significantly affect structural integrity 


and electrochemical performance of the battery. 


2.2 Theoretical and Numerical Framework 


We describe the multi-particle behavior of a composite electrode by considering various 
scales. We consider first a representative volume element (RVE) of the electrode mi- 
crostructure through a continuum description. Subsequently, we focus on the behavior of 
individual active particles and finally transition to the behavior of interfaces. This descrip- 
tion is shown schematically in Fig. 2.1 with a focus on the treatment of electrochemical 


behavior of the composite electrode. 
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Figure 2.1: Schematic description of the theoretical and numerical framework employed 
with a focus on electro-chemistry. (a) Continuum description of the microstructure RVE 
with N particles and galvanostatic conditions enforced. (b) Particle description showing 
discretization of current over M elements over each particle surface. (c) Description of the 
surface element in which non-linear reaction kinetics relate local chemical potential jumps 
at each node pair to the current density at each element. Reproduced with permission from 
[17]. 


2.2.1 Continuum Microstructural RVE Description 


We consider first the continuum behavior of a representative volume element (RVE) of the 
electrode microstructure as shown in Fig. 2.1(a). A critical aspect of capturing the electro- 
chemo-mechanical behavior of the electrode microstructure is proper implementation of 
galvanostatic charging conditions. Here, total current / over the entire domain is held 
constant, while current density over a particular particle, ¿(*), is free to fluctuate across the 


particle surface. To ensure galvanostatic charging conditions, we must enforce 
N 
r=) 10, with IU = | i®dA (2.2.1) 
A Alk) 


where N is the number of particles in the RVE. Eqn. (2.2.1) serves as a constraint, which 
must be numerically enforced over the entire simulation domain. We assume from the 


onset that voltage across all particles surface is constant. We thus restrict ourselves to 
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physical problems in which transport of ionic species through the matrix is sufficiently fast 
and variations in electric potential through the electrolyte are negligible. Consistent with 
this assumption, electric potential across the microstructure RVE is then constant through- 
out. From a numerical perspective, this has the added benefit that we may exclusively 
focus on the interfacial kinetics and active particle behavior without having to resolve ion 
migration and the associated electric field across the conducting matrix. We note that the 
proposed numerical framework may be extended to remove the constraint that voltage be a 
constant across all particle surfaces by considering the physics of ion migration through the 
RVE conducting matrix. The development presented here is a simplification of this more 
general case, but should serve equally useful in the development of future models, which 
incorporate this additional physics. 

To model non-equilibrium interfacial reaction kinetics, we employ the phenomenolog- 
ical Butler-Volmer equation. At every point on the surface of a particle, current density, 2 


is related to the overpotential, 7 through (cf. Bazant [126]) 


¿ato (e (SED) e (252), a 


Here, ¿p(c) is the concentration dependent exchange current density, R - the gas constant, 
Y - the absolute temperature, F - the Faraday constant, and a - a symmetry factor, repre- 
sentative of the fraction of surface overpotential promoting anodic or cathodic reaction at 
the electrode interface, such that 0 < a < 1. The exchange current density, ¿p(c) depends 


on the species concentration through 
io(@) = Fko(1 ae, (2.2.3) 


where ko denotes the reaction constant, while ¢ represents the normalized species con- 
centration detailed in Sect. 2.2.3. From thermodynamics, we may relate the overpoten- 


tial, y to a jump in chemical potential across the particle-electrolyte interface through 
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n = Ap — Af = —Ap/F, with Ay the jump in chemical potential of the diffusing 
species at the reaction site of interest *. This definition of 7, along with the specialization 
a = 0.5, which is typical of elementary single-electron transfer reactions, allows us to 
write (2.2.2) as 
i = 2i9(C) sinh Ga : (2.2.4) 
Importantly, this yields a relationship for the local current density at a point on the particle’s 
surface as a function of the chemical potential jump at the same location. Finally, for the 
microstructure RVE as shown schematically in Fig. 2.1(a), we prescribe periodic boundary 
conditions. 
Particle description and numerical discretization: Consider now a single particle 
discretized with finite elements as shown schematically in Fig. 2.1(b). The behavior of this 
particle can be described through its interfacial and bulk behavior. For clarity, we discuss 


these behaviors separately. 


2.2.2 Chemo-Mechanical Interfacial Particle Behavior 


Focusing first on the chemo-mechanical interfacial behavior, the particle is discretized with 
finite elements resulting in M elements on the surface of the particle. Galvanostatic condi- 


tions are now enforced by writing the constraint (2.2.1) in discretized form as 


T= DA with J) = Nene (2.2.5) 

k=1 p=1 
Here, 1(*) is the total current over the k-th particle, while i”) and a(*») denote the current 
density and area of the p-th surface element of the k-th particle. It is important to recall that 
we allow for and expect current density variations throughout the surface of the particles. 
Hence, we make no assumptions about current density at any point of the particles surface, 


rather we only enforce the galvanostatic constraint that the total current, J be a constant. 


‘Consistent with [42], we define y = (1/F) (uzi — (rit + He- )) 
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To enforce the galvanostatic constraint — and motivated by the physical nature of the 
problem containing non-linear reaction kinetics at the particle/matrix interface — we in- 
troduce surface elements to discretize this interface as schematically shown in Fig. 2.1(c). 
These elements (analogous to traditional cohesive elements in fracture mechanics) are zero- 
thickness elements, which can be thought of as lines in two-dimensional simulations or 
surfaces in three-dimensional simulations. Note that the nodes (and hence degrees of free- 
dom) of the bulk continuum elements used to discretize the particle and the nodes on one 
side of the surface elements are shared. 

We enforce the constraint that voltage across all particles be a constant by constraining 
the chemical potential degree of freedom on the surface element adjacent to the matrix to 
be equal to —V F, with V df (V — Vo), V the electrode voltage, and Vo a reference voltage. 
That is, in Fig. 2.1(c), if we envision that u and uz are connected to the active particles, 
and u3 and u4 are connected to the conducting electrolyte matrix, we restrict our surface 
elements such that u3 = u4 = —V F (c.f. Di Leo et al. [42]). 

We may now write the reaction kinetics (2.2.4) for the current density, i(*P) at an ele- 


ment on the surface in a discretized form 


1 Au EP) 
i) — 2 sinh G a (2.2.6) 


with Ap) - the chemical potential jump for the p-th surface element on the k-th particle. 
Note that the concentration-dependent exchange current density, ier (E) is computed at 
each surface element in order to account for variations in concentration along the particle 
surface. The chemical potential jump, Aj”) is computed using numerical interpolation 
between the chemical potential jumps Au = u4 — uy and Ay?) = uz — u2, details of 
which can be found in Appendix A.2. 


The mechanical behavior of the surface elements is also critical as it allows one to 


capture the manner in which the particle-matrix interface can delaminate due to mechanical 
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stresses. The mechanical behavior relates displacement jumps, Au, across the interface to 
the associated tractions t — in a fashion analogous to how jump in chemical potential, Ay 
relates to the current density, 2 across the interface. Broadly, we allow for a finite (yet large) 
stiffness of the interface, which relates the separation of the surface element nodes to the 
tractions (stresses) generated at the interface. Damage initiation and evolution criteria are 
then prescribed to capture interface delamination. 

Restricting ourselves to a two-dimensional formulation, the normal ty, and tangential 


tr, tractions at the interface prior to damage initiation are computed through, 
tn = KyAuy, and tr = KrAur, (2.2.7) 


with Auy and Aur denoting the normal and tangential separations respectively. Note that 
in Fig. 2.1(c), for brevity, we have only schematically depicted the normal displacement 
jump (i.e. separation). Following Camanho and Da Vila [127], damage initiation at the 


interface is modeled following a quadratic traction function of the form, 


2 2 
(22) +(2) =1, (2.2.8) 
N T 


where tẹ and tF represent the normal and tangential cohesive strengths (i.e. the peak value 


these tractions may attain when the deformation is purely normal or purely tangential to 
the interface), and (-) denotes the Macaulay bracket. 

To model the evolution in mechanical damage, a scalar variable, Dmech 18 introduced, 
which denotes the overall state of mechanical damage for an element on the surface of the 


active particle. As the state of damage at the interface evolves, tractions decay linearly with 
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mechanical damage and are computed through 


(1 = Dmecn) K nAun, if KyAun > 0, 
ty = 
KyAuy, otherwise (to avoid surface interpenetration), (2.2.9) 


tr = (1 = D mech) KpAur. 


To model the evolution of damage, it is common to introduce an effective displacement 
jump, dm = ((Auy)? + Au?.)'/?, which concurrently accounts for both normal and tan- 
gential opening at the interface. We may then model the linear softening of tractions due to 


mechanical damage through an evolution equation for Dec of the form 


On — 0) 
OR (Om — Jn) 


mech — 


Doae lOs (2.2.10) 


Here, 6°, is the effective displacement jump at damage initiation, ôf, is the effective dis- 
placement jump at complete failure, and 6™** is the maximum value of the effective dis- 
placement attained during the loading history, a monotonically increasing quantity. Crit- 
ically, this enforces that damage processes are irreversible in nature, namely once the 
particle-electrolyte interface delaminates, it can not heal. 

Considering positive normal displacement jumps, the effective displacement jump at 


damage initiation, 6°, may be computed as 


1 + 8? pr tE qu 2 pr 2 a 
0 _ 50 50 — NT | 2 T | 2 N 
ll a Ky KN or? (GE) + (aE) 


In writing (2.2.11), we have made use of the fact that 9, = t%/Ky and 6%. = t¥/Kr. 


Additionally, 8 = Aur/Auy denotes the mode-mixity parameter, which accounts for the 
different contributions sourcing from mode I (normal) and mode II (tangential) loading 


of the active particle-electrolyte interface. Similarly, the effective displacement jump at 
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complete failure, ôf, may be computed as 


2 2 =1 
a ) (37) | (5) (2.2.12) 


with Gy and Gr, the normal and tangential fracture energies. In writing (2.2.12), we 


further specialize our model such that K = Ky = Kr. 

As shown from (2.2.11) and (2.2.12), the material properties for this model are the 
stiffnesses, presumed to obey K = Ky = Kry, the fracture energies {G y, Gr}, and the 
critical tractions for damage initiation {t{, 17). Of these parameters, K may be prescribed 
arbitrarily provided it is large enough not to introduce artificial compliance in the simu- 
lation. The physical material parameters are thus only the critical normal and tangential 
tractions {t{\,¢#}, and the fracture energies {Gr,Gy}. Note that the mode-mixity pa- 
rameter, 8 = Aur /Auy depends on the loading path and is not a material property. The 
limiting case 6 — 0 corresponds to pure normal separation, while the limiting case P — oo 
corresponds to pure shear. 

To impart stability to our numerical solver, we employ a viscous regularization scheme 
to model the evolution of mechanical damage at the interface. Following Hamitouche et 
al. [128], we introduce a small viscosity, \ (small with respect to the characteristic time 
increment) and compute the viscous mechanical damage as 

mech = Dmech — va (2.2.13) 


mech — ô 
m 


Henceforth, when referring to mechanical damage, we are employing the viscous damage 
definition above. Viscosity values employed in all simulations in this work ranged between 
0.1 and 0.3s. 


Finally, to model the effect of mechanical damage on reaction kinetics, current density 
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across the interface, viz. (2.2.6), is penalized through 


1 Ape) 
jhe) — (1 — Dat” sinh G a 


J with Denem € [0, 1], (2.2.14) 


where Denem denotes the chemical damage. We specialize the relationship between Dinech 
and Denem in Sect. 2.4.2 and Sect. 3.1 when we specialize our models to a particular phys- 
ical system. As captured through (2.2.14), this model allows one to capture the loss of 
chemical reactivity at a point on the surface of an active particle due to concurring mechan- 
ical damage. The effect of this coupling will be studied in detail in the simulation results 


shown in this thesis. 


2.2.3 Bulk Particle Behavior 


We present here, for completeness, a brief summary of the coupled chemo-mechanical, 
diffusion-deformation theory for modeling the bulk behavior of active particles. Details of 
the theoretical framework may be found in Di Leo et al. [42]. This framework models cou- 
pling of species diffusion with large elastic-plastic deformations due to volumetric changes 
caused by the diffusing species. The active particle is treated as a homogeneous mixture 
of the active material and Li, with the molar concentration of Li per unit reference volume 
denoted by cz and per unit current volume by c = J7*c,, with J = det F, and F the de- 
formation gradient. We further define a normalized concentration € = Cg/Cmax € [0, 1]. We 
employ the decomposition F = F*F”F* of the deformation gradient into elastic, plastic 


and swelling parts, where 
F*=(J9% 41, with J°=1+0(q— Co), (2.2.15) 


is the swelling distortion, (2 is a constant partial molar volume of the intercalating Li in the 


host material, and cy is the initial concentration of Li. Further, we assume J? = det F? = 1, 
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such that plastic deformation be volume conserving. The Cauchy stress is given by 


2 


PEA (aor + (xo E 


co) (wEs)1) ; (2.2.16) 
where G(c) and K (c) represent the concentration-dependent shear and bulk moduli respec- 
tively, while Ey, denotes the spatial logarithmic elastic strain. To complete the mechanical 
portion of the theory, plastic distortion is taken to evolve according to 


; o [3M8 
F? = DPF", wih D?=-@ (55) , e > 0, F?(X,0) = 1, (2.2.17) 
o 


_ def 
where = ,/3/2|M5| defines an equivalent tensile stress, and € denotes an equivalent 


tensile plastic strain rate. 
Turning attention to diffusion, spatial flux j, of the intercalating Li is taken to depend 


on the spatial gradient grad u, of the chemical potential through 
Do - n 
j=-—mgeradu, with m= pg —C)>0__ the mobility. (2.2.18) 


Here, Do is a constant diffusion coefficient. The chemical potential of Li in the active 
particle is given by, 


w=p°+RIIn (75) — OM, (2.2.19) 
= E 


where we have have defined an activity coefficient y through RV In(y) = Xipan- n- 


¿m—1) 


, whose coefficients a, are fitted to experimental open circuit potential data. Im- 
portantly, the chemical potential, (2.2.19) is stress-coupled through the Mandel stress, M° 
which may be related to the traditional Cauchy stress, T through an elastic rotation and 
elastic volumetric scaling [42]. 

The above fields are governed by two partial differential equations expressed in the de- 


formed body. As is standard, the mechanical problem (displacement field) is governed by 


a local force balance, divT = 0, with appropriate boundary conditions requiring that dis- 
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placement or traction be prescribed at a given point on the exterior of the body. The lithium 
diffusion problem (chemical potential field) is governed by mass balance, ¢, = —J divj, 
where appropriate boundary conditions require one to either prescribe the chemical po- 
tential or the spatial flux at a given point on the surface of the body. Finally, we note that 
current density, 2 as described in Sect. 2.2.2 and normal flux of species 7 = j-n, with n, the 
normal to the surface of an active particle at a point, are related simply through j = —¿/F. 
Details of the numerical implementation, with a focus on the implementation of the chemo- 
mechanical surface elements, are presented in Appendix A following the work of Di Leo 
[129] and Chester et al. [130]. 

Next, in Sect. 2.3 we illustrate the capability of the proposed framework to (1) model 
galvanostatic conditions and (11) capture non-uniform current distribution across particles 
along with the role of mechanical delamination at the particle surface through a simple 


demonstration of a two-particle composite electrode RVE. 


2.3 Numerical Framework Demonstration 


To illustrate the capabilities of our numerical framework, we present here a simple sim- 
ulation of a two particle system. Both particles are embedded in an elastic electrolyte 
matrix and their surfaces are discretized with the chemo-mechanical surface elements de- 
tailed in Sect. 2.2.2 above. As described by (2.2.14), we allow here for mechanical degra- 
dation of the interface due to normal stresses leading to a loss in local reactivity. The 
particles are modeled as chemo-elastic (no plastic deformation) and incur 2% volume de- 
crease when fully lithiated. The specific material parameters are the same as those used 
in Sect. 3.1 for a LCO-LGPS composite electrode system. For this illustration, we im- 
pose galvanostatic charging conditions at a C-Rate of 1C, with the additional aforemen- 
tioned constraint that voltage over the particles surface be a constant. Mechanically, we 
impose periodic boundary conditions on the electrode domain, effectively treating it as 


a representative volume element (RVE). For illustration purposes, we induce significant 
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variations in current distribution between particles by assigning the left particle a reac- 
tion constant, ky = 9 - 1076 mol/(m?s), while the right particle has a reaction constant, 
ko = 9 - 1077 mol/(m?s). 

Fig. 2.2(a) shows the total current flowing to the left (blue line) and right particle (or- 
ange line), as well as into the entire simulation domain (yellow line) as a function of time. 
We observe that we can successfully maintain galvanostatic charging conditions, while cur- 
rent distribution to the two particles varies in time. Until time t = tą, no interfacial damage 
occurs and the total current is distributed such that more current (i.e. charge) flows to the 


left particle owing to its faster reaction kinetics. 
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Figure 2.2: Example of a two-particle numerical simulation. (a) Total current distribution 
to both particles as well as the entire simulation domain, illustrating the ability to cap- 
ture galvanostatic charging conditions. (b - c) Local current densities as a function of the 
normalized distance (s1, S2) along the particle surface taken at time t=t, and f=f,. (d - e) 
Contours of maximum in-plane principal stress, cı over the matrix and normalized concen- 
tration, ¢ over the particle domain taken at time time t=t, and t=te. Also illustrated in white 
are regions of decohered interface. Reproduced with permission from [17]. 


After t=ta, interfacial damage occurs over both particles, with mechanical decohesion 
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starting in the left particle first owing to its higher extent of lithiation (and hence higher 
mechanical stresses at the interface). Mechanical degradation now causes current to redis- 
tribute over both particles, as dictated by the interplay of two mechanisms: 1) the variation 
in reaction constant among particles, and 11) the state of mechanical damage on the particle 
surface. As interfacial damage across both particles evolves while ta < t < tp, the inter- 
play of these mechanisms first leads to an increase in current flowing to the left particle, 
while current flowing to the right particle decreases to maintain the galvanostatic charging 
condition. 

Fig. 2.2(b) shows the distribution in current density over each particle as a function 
of the normalized contour distance (i.e perimeter) of the particle at time t=t,. The start 
of the normalized contour variables is noted in Fig. 2.2(d) as sı and s2. We can clearly 
observe the inherent coupling between mechanical degradation and interface kinetics (1.e 
current density). Here, we observe regions over the particles surface with zero (or near zero) 
current density. These regions, marked by white arrows in Fig. 2.2(d), represent elements 
at the surface of both particles which have completely decohered and are now chemically 
isolated and unable to sustain reaction. In Fig. 2.2(d), we show contours of normalized 
concentration, C, over the particle simulation domains, and contours of the maximum in- 
plane principal stress, dv, over the electrolyte matrix at time t = tp. 

After time t=f,, a significant portion of the surface of the left particle has decohered and 
can no longer sustain sufficiently large currents. Thus, there is a decrease in current to the 
left particle as shown in Fig. 2.2(a). In turn, to preserve a state of galvanostatic charging, 
the right particle sees an increase in total current. Fig. 2.2(c) illustrates the distribution in 
current density over each particle as a function of the normalized contour distance at time 
t = t,. Across both particles, we observe that most of the interface can no longer carry 
any significant current (i.e. ability to sustain reaction). In turn, the interface portions which 
remain intact carry a higher current density when compared to Fig. 2.2(b), since the same 


amount of total current prescribed to the electrode domain needs to now distribute over a 
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smaller active area available for reaction. Fig. 2.2(e) shows the corresponding simulation 
results at t = t., where we can observe significant damage accumulated now over the 
surface of both particles. 

This simple numerical example serves to illustrate the ability of our numerical frame- 
work — enabled by the development of chemo-mechanical surface elements — to cap- 
ture: 1) the chemo-mechanical interactions between particles under galvanostatic charging 
and the associated non-uniform current distribution to different particles, and ii) the non- 
uniform variation in local current density which develops across each particle. 

In Sect. 2.4 to follow, we deploy this numerical framework to study mechanical integrity 
and electro-chemical performance in composite electrode systems of relevance to liquid 


LIBs. 


2.4 Multi-particle interactions in double walled a-Si nanotube electrodes 


We consider now a multi-particle electrode composed of hollow double-walled a-Si nan- 
otubes, which are chemically connected, but mechanically isolated, representative of the 
limiting case of a liquid LIB in which mechanical stress transfer through the electrolyte 
is negligible. The simulation mimics the experimental research of Wu et al. [116], who 
manufactured a new anode architecture composed of an ensemble of double-walled a-Si 
nanotubes with the particular objective of preventing failure of the solid electrolyte inter- 
phase (SEI), see Fig. 2.3(a). As has been studied in the literature, volumetric expansions of 
active particles can induce large mechanical stresses on the SEI layer, causing it to fracture 
and fail [44, 131, 132]. The double-walled a-Si nanotubes proposed by Wu et al. [116] 
consist of an a-Si hollow nanotube, whose exterior is oxidized to form a mechanically stiff 
SiO, layer, thus forming a double-walled structure. Due to the high relative mechanical 
stiffness of the exterior SiOz layer when compared to the softer a-Si core, volumetric ex- 
pansions are accommodated primarily through expansion into the hollow internal cavity, 


and minimal strains and stresses are incurred by the SEI layer forming on the exterior 
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surface. In our previous research efforts [42], we investigated the chemo-mechanical be- 
havior of a single hollow double-walled nanotube to gain insight on the coupling between 
mechanics and chemistry in such a confined geometry. However, we did not address the 


multi-particle nature of the electrode or the potential for mechanical damage. 
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Figure 2.3: Multi-particle modeling of hollow double-walled nanotubes. (a) Shows a rep- 
resentative SEM image from the experiment of Wu et al. Reprinted with permission from 
[116]. (b) Shows a schematic three-dimensional representation of a single tube denoting 
the simulation domain as a single sliver of the tube in the axial direction. (c) Shows the 
numerical discretization using finite elements of the a-Si core and the SiO, shell, including 
the discretization of the interface using chemo-mechanical surface elements. Reproduced 
with permission from [17]. 


In this work, each individual tube is modeled through an axisymmetric row of finite 
elements as shown in Figs. 2.3(b) and (c). Here, we assume the length of the nanotubes 
to be much larger than their diameter. This effectively enables one to model each nan- 
otube as a row of elements with the flux of Li-ions occurring in the radial direction only. 
The bottom surface of each row of elements is constrained to have zero displacement in 
the e, direction, while the top surface is constrained to remain flat, but may displace. 
The simulation domain is broken into two parts, namely an a-Si core and the confining 
SiO» shell, with chemo-mechanical surface elements discretizing the interface between the 
two domains, see Fig. 2.3(c). The SiOz shell is prescribed a purely elastic mechanical 


behavior and we neglect transport of Li through this layer. The a-Si obeys the chemo- 
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mechanical framework described in Sect. 2.2.3. Following Di Leo et al. [42], the plastic 
behavior of the a-Si shell is modeled by introducing a positive-valued, stress-dimensioned, 
and concentration-dependent yield strength Y(c) > 0, and assuming a no-flow condition 
of the form a < Y (c). During plastic flow, e” > 0, and the equivalent tensile stress is taken 


to be equal to a rate-dependent flow strength, 


pN 1/m 
AA (=) | with Y(2) = Y + (Yo — You) exp (-2) DA) 


€0 *k 


where Y, > 0 is a positive-valued, stress-dimensioned constant, éy is a reference tensile 
plastic strain rate, and m is a measure of the strain-rate sensitivity of the material. All 
material properties for the a-Si and SiO, are adopted from the work of Di Leo et al. [42] 
and the reader is referred to this reference for numerical values. 

The multi-particle behavior is modeled by including N simulation domains represent- 
ing N particles as shown in Fig. 2.3(c). It is important to note that although these particles 
are mechanically isolated (i.e. that they do not contact each other), they are chemically 
connected through the presence of the chemo-mechanical surface elements. Through these 
elements, we enforce galvanostatic charging conditions across all particles under the as- 
sumption of a constant voltage maintained across the active particle-electrolyte interface. 

The numerical framework described above allows us to investigate the interplay of me- 
chanics and chemistry in multi-particle systems and the associated effect on electrochemi- 
cal performance. First, we investigate the effect of statistical variations in electrochemical 
properties — namely the reaction constant, ko between the active particles and the elec- 
trolyte — on current distribution across particles and interfacial mechanical stress build up. 
Second, we investigate the potential presence and role of interfacial damage on electro- 
chemical performance and capacity fade while cycling. In all of the following studies, we 


consider N = 50 active particles. 
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2.4.1 Electrochemical performance and stress generation with varying reaction rate - ko 


We consider first the effects that would arise in a multi-particle system from having a dis- 
tribution of reaction constants, ko affecting the reaction kinetics at the interface of the par- 
ticles and the electrolyte. Such a variation in kọ may arise experimentally due to a number 
of mechanisms including poor bonding of the a-Si/SiO» interface, failure of the SEI, and 
pre-existing damage of the active material during manufacturing. The reaction constant, 
ko dictates reaction kinetics at the interface (cf. Eq. (2.2.3)), which is modeled through 
the chemo-mechanical surface elements, and is prescribed a log-normal distribution. The 
mean of the distribution is prescribed as the logarithm of the baseline reaction constant, 
ko = 3.25 - 1077 mol/(m?s) reported in Di Leo et al. [42]. The standard deviation (SD) is 
varied in our investigation and we consider values of SD= 1.0 and SD = 2.0. This choice is 
motivated to produce a significant variation in local current distribution over the simulation 
domain in order to understand how these local variations affect the overall electrochemical 
response of the multi-particle system and concurrent stress generation. Consistent with the 
experiments of Wu et al. [116], the simulation domain is cycled under galvanostatic con- 
ditions at a C-Rate of 1C with voltage limits of 0.01 V and 1 V for three half-cycles. We 
note that in this first set of simulation results, we do not allow for mechanical interfacial 
damage, which will be considered in Sect. 2.4.2. 

To understand the effect on electrochemical performance of variations in kọ, we com- 
pute both the local (single nanotube) and global (whole ensemble of nanotubes) voltage 
vs. state-of-charge (SOC) response during cycling. Recall that Voltage is assumed to be 
constant throughout the simulation domain and is an unknown which is solved for in the 
finite element framework in order to ensure the galvanostatic charging constraint is met. In 
generating the Voltage vs. SOC curves subsequently shown, we ignore the results from the 
first half-cycle in our simulations. 

The simulation results are shown in Fig. 2.4. First, we focus on the individual behavior 


of three of the 50 nanotube particles from the SD=1.0 simulation, which have been chosen 
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to represent the lowest, average, and highest reaction constants present in the simulation, 
and differ in ko value by approximately one order of magnitude. An important feature of 
multi-particle simulations evident in Fig. 2.4(a) is the “smoothing” of the Voltage vs. SOC 
behavior of the individual particles as cycling is reversed. This feature is enhanced as the 


reaction constant, ko is decreased, as shown in Fig. 2.4(a). 
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Figure 2.4: Simulation results for a 50-particle system with a statistical distribution in re- 
action constant ko. (a) Local (individual nanotube) voltage vs. state-of-charge behavior for 
three particles spanning the range of ko values considered. (b) Current vs. time response 
for the three particles under consideration compared against the total current in the simula- 
tion domain. (a) and (b) are for the case of SD = 1.0. (c) Global (entire simulation domain) 
voltage vs. state-of-charge behavior for a single particle, simulations with ko distributions 
with SD = 1.0 and SD = 2.0, and the experimentally measured response of Wu et al. [116]. 
Reproduced with permission from [17]. 


This behavior is tied to chemical interactions between particles with varying ko, which 
in turn impact local current distribution to individual particles as shown in Fig. 2.4(b), 
where we compare the current applied to the entire simulation domain (purple line, right 
y-axis) with the local current distributed respectively to the three nanotube particles under 
investigation (left y-axis). While the entire simulation domain experiences a sharp change 
in total current during cycle reversals (purple line), the individual particles experience a 
gradual variation in current, which depends on their local reaction kinetics. Particles with 
sluggish (low) reaction kinetics experience a gradual reversal in current, which in turn 
affects the contribution of the overpotential to the total Voltage vs. SOC behavior and 


results in “smoothing” of the Voltage vs. SOC behavior as shown in Fig. 2.4(a). 
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The global behavior of the system is shown in Fig. 2.4(c). Results are shown for 
two simulations with log-normal ko distributions with a standard deviation of SD=1.0 and 
SD=2.0, as well as for a single-particle (standard deviation of zero). Also shown is the 
experimentally measured response of Wu et al. [116]. The “smooth” transition in volt- 
age during cycle reversals can be observed in the experimental results shown and is also 
captured when a significant variation in kọ is introduced in the multi-particle simulations 
developed in this work. As shown in Fig. 2.4(c), this feature cannot be captured when 
modeling a single particle system, as the jump in current during cycle reversal will always 
lead to a sharp jump in voltage. This illustrates one of the benefits of developing simulation 
tools which accurately model galvanostatic charging conditions in multi-particle electrode 
systems. Further, variation in local current distribution (see Fig. 2.4(b)) also has a sig- 
nificant impact on the generation of stresses and potential for initiation and evolution of 
mechanical damage as investigated next. 

To investigate the interplay between chemical performance and mechanics, we focus on 
the generation of normal stresses at the a-Si/SiO2 interface. Here, we demonstrate the man- 
ner in which variations in local reaction kinetics of each nanotube — as prescribed through 
a variation in the reaction constant ky — affect the generation of mechanical stresses at 
these interfaces. Fig. 2.5(a) shows the normal (i.e. radial direction in Fig. 2.3) component 
of stress o, = Ty at the a-S1/SIO, interface as a function of the local state-of-charge of each 
particle for the same three particles discussed before. As can be observed, the interfaces 
develop different levels of normal stress, which are largely correlated to the local state-of- 
charge of each particle. The local SOC of each particle in turn differs due to variations in 
reaction constant, which promote a non-uniform state of lithiation across particles during 
cycling. 

This behavior is a direct outcome of the elastic-plastic constitutive behavior used to 
model the a-Si core and the large deformation kinematics, which result in an increase (ei- 


ther negative or positive) of the normal interfacial stresses as the a-Si core becomes further 
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Figure 2.5: Interfacial normal stress, Ty at the a-Si/SiO. interface as a function of local 
(nanotube) SOC. (a) Results for three representative particles spanning the range of ko 
values considered at a C-Rate of 1C. (b-c) Results for a particle with a reaction constant 
of ko = 3.25 - 1077 mol/(m?s) at three different C-Rates, with (b) mechanical strain-rate 
sensitivity taken into consideration, and (c) mechanical strain-rate sensitivity neglected. 
Reproduced with permission from [17]. 


lithiated. Importantly, differences between these curves at a given state-of-charge arise 
from the mechanical strain-rate sensitivity of the material. Since particles experience dif- 
ferent currents, they also experience different loading rates, which in turn changes their me- 
chanical response. In summary, the magnitude of normal interfacial stresses is controlled 
both by the total amount of Lithium reacted into the a-Si core and the rate of lithiation. 
From an interfacial mechanical damage perspective (which will be further investi- 
gated in the next section), this poses an interesting question, which is to investigate the 
C-Rate sensitivity of interfacial mechanical stresses. Fig. 2.5(b) shows the normal inter- 
facial stresses as a function of the local SOC for the nanotube with a reaction constant 
ko = 3.25 - 1077 mol/(m?s). Three simulations are shown with C-Rates of 1C, 1/10C, and 
1/20C. We can observe the counter-intuitive result that a decrease in C-Rate does not nec- 
essarily lead to a decrease in the maximum (positive or negative) normal interfacial stresses 
experienced at the a-Si/SiO, interface. The rather constant nature of the maximum normal 
interfacial stresses during cycling shown in Fig. 2.5(b) is due to both the extent of lithiation 
and the mechanical strain-rate sensitivity of the material. At higher C-Rates, we achieve 
a lower local SOC for each particle, which would in turn yield lower interfacial stresses. 


However, mechanical strain-rate sensitivity is amplified at higher C-Rates, which in turn in- 
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creases interfacial stresses. These two competing mechanisms for the particular case of an 
a-S1/S102 double-walled nanotube counteract each other to yield a fairly constant normal 
interfacial stress with varying C-Rates. 

To further illustrate this competition, Fig. 2.5(c) shows the same result where the con- 
stitutive behavior of the a-Si core has been modified to have no strain-rate sensitivity by 
setting the parameter Y, = 0 in (2.4.1). As shown in Fig. 2.5(c), normal interfacial stresses 
are now entirely governed by the local state-of-charge of the particular particle and we can 
observe that lower C-Rates always lead to higher interfacial stresses. 

It is important to note that this behavior — which is counter-intuitive to the usual ex- 
perience that higher C-Rates are always mechanically detrimental — arises from the fact 
that here, we focus on the interfacial mechanical stresses at the a-Si/SiO» interface, rather 
then the bulk stresses developed in either the core or the shell. Bulk stresses generally 
arise due to large gradients in concentration or volumetric expansion under mechanical 
constraint. However, owing to the nano-metered sized geometry of this particular anode 
design, concentration gradients are negligible and bulk stresses in this particular geometry 
arise from the mechanical constraint imposed by the SiO, shell on the a-Si core, which is 
no longer free to expand. An in-depth discussion on the single particle mechanics of this 
anode geometry can be found in Di Leo et al. [42]. 

The results shown in this section illustrate the complex nature of multi-particle behavior 
when taking into consideration potential variations in local reaction kinetics by varying 
the reaction constant, ko through a statistical distribution. The results demonstrate that 
variation in kọ can significantly affect local current distribution among the various active 
particles and in turn influence the global voltage vs. SOC behavior. Further, variation in 
ky can also significantly affect the generation of mechanical stresses with a focus on the 
development of normal interfacial stresses at the a-Si/SiO2 interface. These stresses could 
potentially lead to mechanical damage through decohesion of the interface, which would in 


turn affect chemical connectivity of the particles to the electrode. In the following section, 
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we investigate capacity loss in this prototypical liquid LIB electrode due to mechanical 


damage and loss of connectivity through use of the developed surface elements. 


2.4.2 Capacity loss due to mechanical interfacial decohesion 


We now shift our attention to modeling capacity loss due to mechanical interfacial de- 
cohesion in the electrode under consideration. Unlike the previous section, we consider 
now a statistical variation in mechanical properties associated with the cohesive strength 
of the a-Si/SiO2 interface, while maintaining the reaction rate constant at kọ = 3.25 - 
1077 mol/(m?s). The choice of varying the cohesive strength +5, as opposed to the reaction 
rate ko, in this section is made for pragmatic reasons. The material property, ko has been 
previously experimentally determined in [42], while the specific value of tẹ, as will be dis- 
cussed below, is unclear from the literature and can vary over a significant range depending 
on manufacturing conditions. 

Across the simulation domain containing N = 50 active particles, we prescribe the nor- 
mal cohesive strength t% using a normal distribution with a mean strength of 160 MPa and 
a standard deviation of 70 MPa. The mean cohesive strength is chosen to approximately 
match with the average normal stresses at the interface across an ensemble of nanotubes 
with varying ko as determined in Sect. 2.4.1. The standard deviation is such that at least 
50% of the particles will incur mechanical damage during cycling at different C-Rates. We 
consider only normal mechanical damage (i.e. 6 = 0 in Eqns. (2.2.11) and (2.2.12)) with 
ty as defined above, and simply take tọ in (2.2.8) to be a large number to guarantee no 
damage initiation due to tangential interfacial stresses. The prescribed variation in cohe- 
sive strength is intended to showcase the utility of our modeling framework to study the 
interplay of chemo-mechanical interactions and role of damage on electrochemical perfor- 
mance and capacity fade through consecutive cycling. The normal stiffness of the interface 
is taken equal to Ky = 6 - 10% MPa/m, which is sufficiently large not to introduce artifi- 


cial compliance in the simulation domain. To define damage evolution, we use a baseline 
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fracture energy of Gy = 0.25 J/m?. In combination, these material properties define all the 
necessary parameters for the constitutive behavior summarized in eqns. (2.2.8) to (2.2.10). 
Finally, for simulations in this section, we consider chemical damage to evolve equal to 
mechanical damage, that is Dehem = Dmech- 

The simulation domain is cycled under galvanostatic conditions with varying C-Rates 
and voltage limits of 0.01 and 1 Volts. We simulate 27 cycles of the electrode with a 
baseline C-Rate of 1/2C and alternating cycles at higher and lower rates. The complete 


sequence of C-Rates is given by, 


C-Rate € [1/2C, 2C, 1/2C, 4C, 1/2C, 1/4C, 1/2C, 1/10C, 1/2C] (2.4.2) 


where three cycles are spent at each particular C-Rate. This sequence allows us to asses 
how damage incurred at higher or lower C-rates — with respect to the baseline 1/2C rate 
— affects electrochemical performance. By computing the capacity of the electrode at a 
C-Rate of 1/2C before and after cycling at a different rate, we can understand how damage 
incurred at a particular rate affects the capacity of the electrode. Following Di Leo et al. 
[42], a maximum capacity of 3.579 Ah/g for Si is used to convert SOC data to an equivalent 
capacity. 

Simulation results are shown in Fig. 2.6. The top row shows the voltage vs. SOC 
response for the entire simulation domain, while the bottom row shows the simulation do- 
main, where in the a-Si domain we show contours of the normalized concentration, c. The 
SiO» shell is shown as a grey region without contours. Fig. 2.6(a) shows the electrochemi- 
cal response of the electrode at the end of the third conditioning cycle at a C-Rate of 1/2C at 
which point the battery accrues a significant amount of damage. Fig. 2.6(c) highlights with 
red arrows the particles which have completely decohered from their SiO shells during the 
conditioning cycles and are now chemically isolated from all other particles. 


Due to mechanical decohesion, the electrode stabilizes to a new electrochemical win- 
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Figure 2.6: Simulation result of cycling a multi-particle electrode composed of double- 
walled, hollow nanotubes. Figs. (a) and (b) show the Voltage vs. SOC response of the 
electrode where (a) compares a simulation with and without damage after the first three 
conditioning cycles and (b) shows results of the simulation with damage after 9, 12, and 
15 cycles. Figs. (c) and (d) show contours of normalized concentration in the a-Si core 
at the end of the 3rd an 15th cycle. The red arrows highlight those particles which have 
completely decohered from their shells and are now chemically isolated. Reproduced with 
permission from [17]. 


dow with a significantly reduced capacity (blue line) compared against a simulation without 
damage (yellow line). Fig. 2.6(b) shows the voltage vs. SOC response of the a-Si anode at 
the end of the 9th (blue line), 12th (orange line) and 15th (purple) cycle during which the 
electrode is intermittently charged at a C-Rate of 4C. By comparing the voltage vs. SOC 
behavior during cycling at the baseline 1/2C rate before (blue line) and after (purple line) 
the intermittent 4C cycles, we can visually see the additional loss in capacity as the electro- 
chemical window shrinks. Fig. 2.6(d), highlights with red arrows the additional nanotubes 
which completely decohere from their SiO» shells during the intermittent 4C cycling. 
Capacity loss due to interfacial decohesion is best characterized through a graph of ca- 
pacity as a function of cycle number as shown in Fig. 2.7. Here, simulation results without 


(blue line) and with (orange line) damage are displayed. The response of the simulated 
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Figure 2.7: Capacity of an ensemble of hollow double-walled nanotubes as a function of 
cycle number during galvanostatic cycling at various C-rates. Results are shown for a 
simulation without interfacial damage (blue line) and a simulation with interfacial damage 
(orange line). Reproduced with permission from [17]. 


electrode without mechanical damage (blue line) in Fig. 2.7 matches the theoretical expec- 
tations for galvanostatic charging within a fixed voltage window, where simulations with 
higher C-Rates have lower capacity and lower C-rates have higher capacity. It is important 
to note that the baseline capacity at a C-rate of 1/2C remains unchanged in the absence of 
mechanical damage and loss of electrochemical contact. The simulated electrode response 
with mechanical damage (orange line) shown in Fig. 2.7 is markedly different. A large 
amount, approximately 75%, of the electrodes capacity is lost in the first three condition- 
ing cycles. This correlates with a large number of nanotubes becoming decohered from 
their SiO» shells as shown in Fig. 2.6(c). During subsequent cycling at 2C, we do not in- 


cur additional damage as we can see that the baseline capacity at 1/2C remains unchanged 
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after cycling at 2C. During cycling at 4C however, the electrode incurs additional damage 
and we can see another significant drop in the baseline 1/2C capacity of the electrode after 
cycling at this rate. This is confirmed by the observations in Fig. 2.6(d), where we see visu- 
ally that additional particles have become decohered. In contrast, no additional damage is 
incurred in the electrode during cycling at C-Rates of 1/4C and 1/10C, which are below the 
1/2C baseline. As shown in Fig. 2.7, the baseline electrode capacity at 1/2C is recovered 
after intermittent cycling at lower C-Rates. 

While use of a different statistical distribution for tẹ would clearly impact these re- 
sults, the qualitative nature of the above study is nevertheless of significant importance and 
demonstrates: 1) the importance of accounting for mechanical damage, which leads to a loss 
of electrochemical activity at the surface of an active particle, and 11) the manner in which 
that can be achieved computationally through use of the developed chemo-mechanical sur- 
face elements. In combination with further experimental calibration of material properties, 
the developed modeling technique should serve useful in analyzing the performance of 
future LIB electrode designs, including chemo-mechanical interactions and role of me- 


chanical damage on electrochemical performance and capacity fade. 
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CHAPTER 3 
ON THE ROLE OF INTERFACIAL MECHANICAL DAMAGE IN COMPOSITE 
CATHODES FOR ALL-SOLID-STATE BATTERIES 


We consider now an application of our modeling framework developed in chapter 2 towards 
composite cathodes for all-solid-state batteries (SSBs), where active particles are connected 
through a relatively stiff solid-state electrolyte (SSE). This chapter focuses specifically on 
a relevant engineering problem where active particles interact both through chemical con- 
nection and direct mechanical load transfer owing to the presence of a stiff solid conductor. 

Common design of composite cathodes incorporates a densely packed architecture of 
active particles, confined on their exterior by a conductive matrix which supports perco- 
lation pathways for ionic transport in addition to complementary additives for electronic 
conduction [16, 24, 133]. The active material undergoes volumetric expansion/contraction 
during insertion/extraction of the ionic species, which, due to the confined nature of the 
all-solid-state architecture can lead to generation of high stresses. Delamination of the in- 
terface between active particles and the SSE constitutes a critical mechanical degradation 
mechanism, which can impact electrochemical performance of composite cathodes and 
SSBs [32, 134]. As a result, efforts aimed at optimizing the performance of composite 
cathodes for application in SSBs require an understanding of interfacial processes, includ- 
ing damage occurring within the electrode as well as the SSE [31, 99, 100]. 

In this chapter, we specialize our framework for a system composed of LiCoO, (LCO) 
active particles and a Lijp)GeP2S12 (LGPS) electrolyte, and integrate our chemo-mechanical 
cohesive element to capture interfacial reaction kinetics and damage. We apply this frame- 
work towards studying the manner in which various material and design parameters affect 
the electrochemical performance of the system. In particular we: i) in Sect. 3.2 investi- 


gate the role of variations in chemo-mechanical properties of the active material and SSE 
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(i.e. volumetric expansion, Young’s modulus), and 11) in Sect. 3.3 investigate the role of 
variations in microstructural descriptors (i.e. particle size distribution, packing density). 
Through this modeling, we elucidate the role of both material properties and microstruc- 
tural composition on the integrity of composite electrodes and overall electrochemical per- 
formance of the solid-state battery. The work presented in this chapter was adapted with 


permission from and published in: 


e Bistri, Donald, and Claudio V. Di Leo. ”Modeling of chemo-mechanical multi- 
particle interactions in composite electrodes for liquid and solid-state Li-ion batter- 


ies.” Journal of The Electrochemical Society 168.3 (2021): 030515. 


3.1 Model Parameters 


For clarity and completeness, we summarize here all material parameters utilized in our 
simulation framework for the LCO-LGPS composite cathode under investigation. Where 
possible, material parameters are taken directly from the literature. The remaining proper- 
ties have been calibrated to experimental data and a summary of the calibration process is 
presented. Material properties for our simulation framework are summarized in Table 3.1 
along with respective sources. In the following subsections, we briefly detail the manner in 


which material properties were chosen or calibrated. 


3.1.1 Mechanical Properties of LCO and LGPS 


The mechanical behavior of LCO active particles is taken as chemo-elastic and modeled 
using the same framework described in Sect. 2.2.3, however here without any plastic de- 
formation. The LGPS electrolyte is treated as a purely elastic solid. Mechanical properties 
of LCO have been investigated in multiple works [135, 136, 137, 138, 139]. However, 
past studies and modeling efforts have resorted to an isotropic assumption, thus neglecting 
variations in mechanical properties of Li CoO% with Li content and crystallographic ori- 


entations. Wu and Zhang [140] applied first-principle calculations in addition to ab-initio 
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tensile tests to study anisotropic and Li concentration-dependent mechanical properties of 
Li,Co0O». Variation in elastic stiffness, shear and bulk modulus with Li content along two 
axes of anisotropy were reported. While our simulation framework does not account for 
anisotropic behavior of Li CoO» in a formal continuum treatment, the aforementioned find- 
ings of Wu and Zhang are incorporated in the model through a Voigt-Reuss-Hill homoge- 
nization scheme, as reported in [141]. Here, a sixth-order polynomial fit is performed to 
capture the variation of Young’s Modulus with Li content from /=108.5 GPa for a pristine, 
unlithiated state to #=252 GPa for a fully-lithiated state. We employ a constant Poisson ra- 
tio of y = 0.22 for LCO. This effectively enables us to capture the experimentally reported 
variation in shear modulus with Li content from G=41.71 GPa for an unlithiated state to 
G=111.38 GPa for a fully lithiated state. A plot of the polynomial fit for both Young’s 
modulus and Shear modulus against experimental data by Wu and Zhang [140] is reported 
in Appendix A.3. 

The mechanical behavior of the LGPS electrolyte is modeled under a linear-elastic 
framework. Several studies have investigated the elastic properties of LGPS through first- 
principle calculations [142, 143]. In Sect. 3.2, we vary the mechanical properties of the 
SSE to investigate the potential impact of changing the SSE composition on electrochem- 
ical performance of the electrode. We use LGPS in our baseline simulations and adopt 
mechanical properties reported by by Wang et al. [142], namely a Young’s modulus of 


E = 37.19 GPa and a Poisson ratio of y = 0.296. 


3.1.2 Chemical Properties of LCO 


Following Di Leo et al. [42], we model the chemo-mechanical behavior of the Li CoO» 
active particles through the framework summarized in Sect. 2.2.3 above. An important 
aspect of using this framework is calibrating the activity coefficient, y through a fitting 
of the open-circuit behavior of the active material. In the absence of mechanical stresses, 


using (2.2.19), one can relate the chemical potential of Li atoms at the electrode surface to 


47 


a stress-free, open-circuit potential U?, according to the following relationship (cf. Bucci 


et al. [115]), 


0 RT E 
i EE ae n(s ) (3.1.1) 


where Uf = —uo/F defines a standard rest potential. Motivated by Verbrugge and Koch 
[144], the concentration dependent activity coefficient y, is given by the following polyno- 
mial representation 

N 


RT In(y) = So ann 2", (3.1.2) 


n=2 

Using eqns. (3.1.1) and (3.1.2), a least-square polynomial regression fit is performed to 
published open-circuit potential data by Mizushima et al. [145] for a Li,CoO2 compound 
against a Li-metal reference/counter electrode to determine the polynomial coefficients in 
(3.1.2). In this work, a seventh order polynomial representation is adopted to determine the 
a, coefficients and Ue with the calibrated values summarized in Table 3.1. A plot of the 
polynomial fit against experimental open-circuit potential data is shown in Appendix A.3. 

To model the diffusion process, a Li diffusion coefficient of Dy = 5.387 - 107% m?/s 
and a maximum molar concentration of Li in the host material Of Cr max = 51555 mol/m? 
are used. These values match with estimates employed in previous modeling efforts for 
a similar electrode material by Wiedemann et al. [146] and Ramadass et al. [147]. We 
note that there is a large variation in reported estimates for the diffusion coefficient of Li 
in Li,CoO» electrodes, with values ranging from 107*% to 10716 m?/s [148, 149, 150]. The 
volumetric changes experienced with added Li content, that is the partial molar volume Q, 
is quantified based on the experimental findings of Reimers and Dahn [151]. It is important 
to note that in contrast to other active materials, Li,CoO2 experiences a contraction upon 
Li insertion, with an approximately 2% contraction in volume reported upon lithiation to 
Lip,gCoO,. Based on these experimental observations, we set the quantity Q = Qemax = 
—2%, which controls the total volumetric change of active particles due to Li intercalation. 


To model interfacial kinetics, we calibrate the reaction constant, ko, to experimental 
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charge-discharge curves for Li CoO» published by Zhang et al. [152]. We do so by ad- 
justing the reaction constant, ko so that energy dissipation during a full cycle (evaluated as 
the area inside the Voltage vs SOC curve) matches with experimental measurements. The 
calibration procedure is repeated at both a C-Rate of C/7 and C/2.7 to ensure consistency 
in calibration at different charging rates. The calibrated reaction constant evaluates to ky = 
9 - 1077 mol/(m?s). A figure of our calibration against the experimental charge-discharge 


curves by Zhang et. al [152] at a C-Rate of C/2.7 is shown in Appendix A.3. 


3.1.3 Interfacial Properties 


To the authors knowledge, no prior data on interface strength has been reported for a 
LCO-LPGS composite cathode and in general experimentally measured interface prop- 
erties within SSBs are scarcely available as pointed out in the review work of Zhang et al. 
[153]. As such, in the numerical investigations in subsequent sections, we prescribe t% 
values over a certain range to provide a qualitative assessment on the effect of damage on 
performance of composite cathodes in SSBs. 

Further, we note that we consider only normal decohesion at the interface (i.e. PB = 0 
in Eqs. (2.2.11) and (2.2.12), and tẹ — oo), since given the limited experimental data 
available, it is not necessary at this stage to include a complex, mixed-mode failure mech- 
anism for the interface. In Sect. 3.2 we vary the cohesive strength between tẹ = 150 MPa 
and ty = 250 MPa. Additionally, in selecting the aforementioned cohesive strengths, an 
effective Young’s Modulus to Cohesive Strength ratio, Esse/t5 € [0,500] is maintained, 
comparable to similar studies in the literature [100]. Fracture energies of Gy = 1J/m? 
and Gy = 1.75J/m? are prescribed for the two different cohesive strengths to main- 
tain a similar normal separation at failure. These values are consistent with the recent 
modeling efforts of Bucci et al. [100]. Stiffness of the interface is prescribed equal to 
K = Ky = 5-10" MPa/m. Choice of a sufficiently high interfacial stiffness is important 


to preserve the correct mechanical behavior of the microstructure prior to damage initia- 
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tion, and eliminate the addition of artificial compliance. The stiffness, K was determined 
by running multiple simulations with increasing values of K and analyzing the develop- 
ment of normal stresses at the interface until a sufficiently large stiffness value to match the 


stress state for a perfectly bonded interface was determined. 


Table 3.1: Material properties for modeling of coupled chemo- mechanical behavior of a 
LiCoO>-Li¡oGeP25:12 composite cathode. 


Parameter Value Source 


Chemical Do 5.387-107% m?/s Wiedemanan et al. [146] 


CR max 51555 mol/m? Ramadass et al. [147] 

Q - CR max -2% Reimers and Dahn [151] 

[a2 — a7|/F [0.1447, 0.4629, -0.7643, Fitted to Mizushima et al. 
2.5326, -3.199, 1.2263] V [145] 

Ug 4.6 V Fitted to Mizushima et al. 


Reaction ko 


9.1077 mol/(m?s) 


[145] 


Fitted to Zhang et al. 


[152] 
Elastic E.ops 37.19 GPa Wang et al. [142] 
ULGPS 0.296 Wang et al. [142] 
Exi,Co0, 108.5 - 252 GPa Wu and Zhang [140] 
VLi,CoO, 0.22 Wu and Zhang [140] 
Interfacial t% [150, 250] MPa This work 
Mechanical Gy [1, 1.75] J/m? This work 
Ky 5 - 101 MPa/m This work 


Finally, as described in Sect. 2.2.2, current density at a point on the surface of a particle 
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is damaged through the parameter Denem, Which is taken to evolve according to 


Donem = 1 — exp(—k - Dmecn) (3.1.3) 


where Dec denotes mechanical damage, which evolves according to (2.2.10), and x = 30. 
This allows us to prescribe a response for mechanical damage which is relatively ductile, 
while the response for chemical damage is of a brittle nature. For the value of k = 30 
chosen, at 15% mechanical damage, Dmech = 0.15, we have already achieved 99% loss of 
chemical conductivity, Denem = 0.99. This eases numerical convergence significantly with- 
out a loss of the ability to capture sudden losses in chemical conductivity due to mechanical 


damage. 


3.2 Effect of varying chemo-mechanical properties on electrochemical performance 


of SSB composite cathodes 


To study the effect of variations in chemo-mechanical properties on electrochemical re- 
sponse of a composite cathode, we consider a two-dimensional microstructure. SEM im- 
ages of a LiCoO2-LijpGeP2S12 composite cathode provided in the work of Zhang et al. 
[16] and shown in Fig. 3.1(a) were used to construct a representative microstructure model 
as shown in Fig. 3.1(b). We model the electrode particles to be elliptical in shape and 
prescribe the distribution in particle major and minor axis to match the experimental SEM 
images. Similar particle dimensions for representative LiCoO, microstructures are also 
reported in the work of Wilson et al. [154]. 

Periodic displacement boundary conditions are prescribed, effectively treating the sim- 
ulated microstructure as a representative volume element (RVE). The surface of each ac- 
tive particle has been discretized with chemo-mechanical surface elements marked by red 
contour lines in Fig. 3.1(b). The particles are electro-chemically coupled by prescrib- 


ing that voltage across the active particles surface remains uniform across as described 
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22 um 


25 um 


(a) (b) 


Figure 3.1: (a) SEM image of a LiCoO>2-Li¡yGeP2S,2 composite cathode microstructure 
from the work of Zhang et al. Reprinted with permission from [16]. (b) Graphical illustra- 
tion of the 2D microstructure mimicking SEM images of Zhang et al. [16] used to model 
the effect of variation in chemo-mechanical properties on electrochemical performance of 
composite cathode. Reproduced with permission from [17]. 


in Sect. 2.2.1. Finally, we note that other complex degradation mechanisms, such as for ex- 
ample the formation of unstable interphases at the active particle/SSE interface can affect 
both mechanical degradation and interfacial kinetics [18]. Treatment of these phenomena 
is beyond the scope of this thesis, where we consider the active particle/SSE interface to 


remain pristine and be affected only by mechanical decohesion. 


3.2.1 Role of SSE stiffness on electrochemical performance and mechanical degradation 


We consider first the role of varying the Young’s modulus of the SSE. We vary the elastic 
modulus between Essg = 20 GPa and Essg = 100 GPa, representing different families 
of inorganic electrolytes, namely sulfides, lithium phosphorous oxynitrides, and garnets 
[31, 153]. The microstructure is cyclically charged-discharged at a C-Rate of 0.5C for five 
half-cycles between voltage caps of 3.8 - 4.6V." 


‘From an experimental standpoint, Li,CoO, is typically charged in the 0.5<x<1 stoichiometric range 
to avoid the detrimental role of phase transitions on battery lifespan. In this work, however, we focus on 
understanding the role of mechanics and microstructural descriptors on electrochemical performance and 
neglect the role of phase transitions on capacity fade. As such, in our simulations, lithiation of Li,CoO2 
particles from a pristine state is conducted. Studies employing a similar stoichiometric range are reported in 
the literature [155, 156, 157]. Additionally, doping techniques for mitigation of phase transitions in LiCoO2 
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First, we establish a baseline behavior by running simulations without interfacial dam- 
age. The results are shown in Fig. 3.2 where (a) shows the average normal interfacial stress 
Ty? experienced in all active particles during charging, while (b) shows the overall Voltage 
vs. SOC behavior for the composite cathode. As shown in Fig. 3.2(a), there is a significant 
increase in interfacial normal stresses (and hence also bulk stresses experienced by the ac- 
tive particles) as we increase the stiffness of the SSE. However, this increase in stress is not 
significant enough to in and of itself cause a significant change in electrochemical behavior 
of the composite cathode as shown in Fig. 3.2(b). This increase in interfacial stresses can 
however have a significant impact on electrochemical performance if one is to account for 


damage and decohesion of the interface as shown next. 
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Figure 3.2: Composite cathode response with no damage and varying SSE Young’s Modu- 
avg 


lus. (a) Average normal interfacial stress, Ty” induced at the active-particle/SSE interface 
vs. SOC. (b) Voltage vs. SOC response of the composite electrode.Reproduced with per- 
mission from [17]. 


The electrochemical Voltage vs. SOC response of the composite electrode with inter- 
facial damage is shown in Fig. 3.3 with (a) tẹ = 250 MPa and (b) tẹ = 150 MPa for the 
baseline case of Q = —2% volumetric contraction. For the case with higher interfacial 
strength in Fig. 3.3(a), we see a decrease in electrochemical window (i.e. capacity) for all 


simulations with a SSE stiffness above 20 GPa. As expected, capacity fade becomes more 


have been successful [158], enabling for cycling at higher voltages. 
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pronounced with increase in SSE stiffness owing to the higher interfacial stress genera- 
tion during deformation of the active particles. As these interfaces become damaged, they 
lose their charge transfer capacity, which has two consequences. First, under galvanostatic 
charging, local current densities must increase to maintain a constant total current, leading 
to an increase in overpotential. Second, some particles may become disconnected from the 
SSE, making them inaccessible for further ion storage. The effect is further exacerbated at 
the lower interfacial strength of tẹ = 150 MPa shown in Fig. 3.3(b), where we note a sig- 
nificant decrease in electrochemical window for all SSE stiffness values considered. These 
observations are consistent with previous modeling efforts by Bucci et al. [100], whereby 


compliant electrolytes are found to perform better at sustaining interfacial integrity. 
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Figure 3.3: Effect of variations in SSE stiffness on electrochemical response of a compos- 
ite cathode with different cohesive strengths ¢{ and volumetric contractions of the active 
particles ©. (a) Simulation results with high cohesive strength ¿$ = 250 MPa and 2% vol- 
umetric contraction, Q = —2%. (b) Results with reduced cohesive strength ty = 150 MPa. 
(c) Results with increased 5% volumetric contraction, Q = —5%. The range of SSE stiff- 
nesses considered are representative of Sulfides, LGPS, LiPON and Garnets. Reproduced 
with permission from [17]. 


3.2.2 Role of active material volumetric change on electrochemical performance and 


mechanical degradation 


We consider next the role of varying the volumetric change incurred by the Li,,CoO, active 
particles during lithiation and its effect on electrochemical performance. Here, we fix the 


normal interfacial strength at tẹ = 250 MPa. This mimics studying the potential effect 
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of changing the composition of active particles to ones with higher capacity and hence 
higher volumetric changes during lithiation/delithiation. We consider a case with increased 
volumetric contraction upon lithiation of 9 = —5%. The effect on the simulated voltage 
vs. SOC response is shown by comparing Figs. 3.3 (a) with Q = —2% and (c) with 
Q = —5%. As observed by contrasting Figs. (a) and (c), increase in volumetric change 
upon lithiation of active particles has a drastic effect on the loss of electrochemical window 
due to mechanical interfacial damage. Specifically, in going from a low to a high stiffness 
SSE, capacity decreases by 56, 84, 90 and 95% when the volumetric contraction is 5%. In 
the extreme case where the SSE stiffness is Essg = 100 GPa and the volumetric contraction 
upon lithiation is Q = —5%, the composite cathode loses almost its entire energy storing 
capacity as the majority of the particles in the electrode undergo complete decohesion from 
the SSE matrix. 

The results in Fig. 3.3 show a homogenized view of the composite cathode by describ- 
ing the voltage vs. SOC behavior of the entire simulation domain. Within the domain, 
there is a complex coupling in the manner in which mechanical stresses develop, damage 
occurs, and interfacial currents re-distribute. This coupling is illustrated in Fig. 3.4 for the 
simulation with Essp = 20 GPa and Q = —5%. The top left image of Fig. 3.4 shows the 
average current density (blue), maximum interfacial normal stress Ty” (orange), and aver- 
age damage (yellow) for “particle 1”, which is labeled in Fig. 3.4(a) as the top-left particle 
in the simulation domain. Figs. 3.4(a)-(d) show two sets of contours. Over the SSE domain 
we show the maximum in-plane principal stress o1, while in the active particles we show 
contours of normalized concentration c. Note that deformations are scaled by a factor of 
10 to better visualize the formation and evolution of interfacial decohesion. 

Fig. 3.4(a) shows a snapshot before any damage has occurred over the simulation do- 
main (corresponding to the first dashed line in the top-left plot). At this point we can 
clearly see the formation of large stresses between particle pairs which are near each other. 


The remaining three contours shown in Figs. 3.4(b) through (d) show the progression of 
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Figure 3.4: Simulation results for the particular case with Essg = 20 GPa and Q = -5%. 
The top-left plot shows the average current density, maximum normal interfacial stress 


Ty”, and average damage for Particle 1 which is shown in (a). Figs. (a) through (d) show 


contours of maximum in-plane principal stress a; in the SSE matrix, while contours of the 
normalized concentration € are shown over the active particles. Note that deformations are 
scaled by a factor of 10. Reproduced with permission from [17]. 


damage as the simulation evolves and the active particles lithiate. Critical to note is the 
behavior of Particle 1 as shown in the top-left corner. At time (b), we observe that Particle 
1 experiences a change in the maximum normal traction it observes as well as the average 
current density over its surface. This is entirely due to crack formation at other particles 
in the simulation domain. This illustrates the complex mechanical and chemical coupling 
which occurs in solid-state composite cathodes. Interfacial damage in one particle, and the 
associated loss in current at the damage site, impacts the mechanical behavior and inter- 
facial current in all other particles in the domain. The phenomena is again illustrated in 
Fig. 3.4(c) where further damage accumulated in the simulation domain — while Particle 1 


remains undamaged — affects the interfacial behavior of Particle 1 as shown on the top-left 
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image. Finally, Fig. 3.4(d) shows the simulation domain at the end of the first half-cycle. 
Here, Particle 1 has incurred significant damage, but remains partially connected to the 
SSE matrix. 

Interfacial mechanical delamination is shown here to play a significant role in the man- 
ner in which the SSE stiffness and choice of active material can affect electrochemical per- 
formance of all-solid-state composite electrodes. Even relatively small variations in these 
material properties can have a significant impact on electrochemical performance. From 
a design perspective, it is thus critical to map out design guidelines of material property 
pairs (such as SSE stiffness and active particle expansion/contraction), which will enable 
for design of solid-state composite cathodes capable of sustaining integrity of the interface. 
This investigation also points out the critical need for further experiments to characterize 
the cohesive strength of different active particle/SSE interfaces as this can be critical to 


modeling the performance of composite electrodes. 


3.3 Modeling the effect of particle size distribution and active material volume frac- 


tion on interface damage in composite cathodes 


We now turn our attention to studying the role of microstructural descriptors (e.g particle 
size distribution, packing density) on electrochemical performance and interfacial integrity 
of composite cathodes in SSBs. These studies provide insight into important design param- 
eters of consideration for composite cathodes and further illustrate the use of the proposed 


numerical framework for modeling the chemo-mechanical behavior of SSBs. 


3.3.1 Role of particle size and size distribution on electrochemical performance and mechanical 


degradation. 


We first study the effect of particle size distribution on electrochemical performance and 
mechanical degradation. As will be shown here, the key factor controlling the role of parti- 


cle size distribution are mechanical particle-particle interactions arising from the presence 
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of the relatively stiff SSE matrix, rather than the actual size (and size distribution) of the 
particles themselves. To illustrate this, we first conduct a simple study of a RVE com- 
posed of a single circular active particle embedded in a SSE matrix with periodic boundary 
conditions as shown in Fig. 3.5(a). Note that we do not allow for interfacial damage in 
these single particle simulations and focus only on investigating the build up of mechanical 
stresses at the interface. Importantly, due to periodic boundary conditions, the simulated 
domain is a repeated RVE and does not represent a single mechanically isolated particle, 
but a repeating unit of particle pairs all of the same geometry. Using the simulation domain 
shown in Fig. 3.5(a), we vary the volume fraction of active material defined through the 
geometric factor f = Rap /Rins, where Rap is the active particle radius and Rins is half the 
length of the simulation domain. The active material volume fraction is varied by chang- 
ing both the particle size and the domain size through i) varying the active particle radius 
Rap € [1, 10] um for a fixed domain size of Rins = 12.5 um, and ii) varying the domain size 


Rins € [6,50] um for a constant active particle radius of Rap = 5 um. 
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Figure 3.5: (a) Single particle representative volume element (RVE) simulation domain. 


(b) Plot of average normal interfacial stress 7," as a function of the geometric factor (i.e 


volume fraction) f = Rap/Rins. Reproduced with permission from [17]. 


avg 


Fig. 3.5(b) shows the average normal interfacial stress Ty * induced in the particle at 
an SOC of 0.8 as a function of the active material volume fraction f = Rap/Rins. Results 
are shown both for simulations varying the active particle radius Rap (circle markers) and 
simulations varying the domain size Rins (x markers). The first critical observation to make 


is that the two simulation results are identical for a fixed active material volume fraction 
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f. This demonstrates that active particle size by itself does not control stress generation at 
the interface, rather it is controlled by the active material volume fraction. In essence, form 
a design standpoint, generation of interfacial stresses is controlled by the electrode active 
material volume fraction (i.e. packing density) - through the presence of particle-particle 
mechanical interactions - and is independent of the actual size of the active particles. Sim- 
ilar results are found in the study by Bucci et al. [100]. 

The second critical observation shown in Fig. 3.5(b) is the importance of particle- 
particle interactions. At low active material volume fractions, roughly below f = 0.3, 
the particle simulated can be considered as “isolated” and we observe that the magnitude 
of INE remains relatively constant and low. As volume fraction of active material increases, 
the particle (through the prescribed periodic boundary conditions) acts in a multi-particle 
behavior and is affected by particle-particle interactions. In this regime of higher volume 
fractions of active material, the magnitude of Ty” increases rapidly with volume fraction 
as particle-particle interactions lead to a continuous increase in interfacial stresses with 
reduction in spacing between particles. 

We now expand our study beyond a single particle to investigate the performance of a 
composite cathode RVE with varying particle size distribution and constant volume fraction 
of active material. Fig. 3.6 shows the simulated RVE domains, all of which have a constant 
active material volume fraction of dam = 30%. Particle size and aspect ratio are seeded 
using a uniform distribution with prescribed lower and upper bounds. We consider active 
particle size bounds in the range of: i) 2 — 6 um, ii) 4 — 5 um, and iii) 6 — 7 um. This allows 
us to investigate the response of the electrode for both a uniform distribution with relatively 
small (4-5 um) and large particles (6-7 um), as well as a distribution with a larger spread in 
particle sizes (2-6 um). 

The RVEs were generated using MicrostructPy [159], an open-source microstructure 
mesh generator, which allows the user to prescribe the active material volume fraction and 


lower and upper bounds for the uniform distribution in particle size as an input. The aspect 
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Figure 3.6: Simulated RVEs with a constant active material volume fraction dam = 30%. 
The particle sizes are given a uniform distribution with lower and upper bounds in the range 
of: i) 2-6 um 11) 4-5 um and iii) 6-74m. The aspect ratio lower and upper bounds fare set 
to 0.9 and 1.5 respectively for all cases. Reproduced with permission from [17]. 


ratio lower and upper bounds for the elliptical particles are set to 0.9 and 1.5 respectively 
for all simulation domains. We cycle the electrode over three half-cycles between 3.8 - 4.6 
V voltage caps at a C-Rate of 1/2C with interfacial damage active. All material properties 
are as described in Sect. 3.1, and for the interfacial cohesive strength we use tẹ = 150 MPa. 

The Voltage vs. SOC response of the three RVEs is shown in Fig. 3.7(a) and compared 
to simulation results with no damage (which are identical for all particle size distributions). 


As shown in Fig. 3.7(a), electrochemical response of the RVE remains consistent across all 
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Figure 3.7: Simulation results for composite cathode RVEs with constant active material 
volume fraction dam = 30% and varying particle size distribution. (a) Voltage vs. SOC 
behavior also compared to simulations with no damage. (b) Evolution of average Denem as 
a function of time for the first half-cycle. (c) Evolution of Voltage as a function of average 
Denem. Reproduced with permission from [17]. 
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RVEs with a fixed dam = 30%, noting again that particle size distribution is not a critical 
factor in determining electrochemical performance. 

The evolution of damage is also consistent across the three microstructures as shown in 
Figs. 3.7(b) and (c). As noted in Sect. 3.1, mechanical damage, Dmecn is taken to evolve 
in a rather ductile manner (so as to avoid numerical convergence issues), while chemical 
damage D nem, defined in (3.1.3), is of a brittle nature. Fig. 3.7(b) shows the evolution 
of average Denem (i.e average accumulated chemical damage) over the entire simulation 
domain as a function of time. We note that all damage occurs during the first charging half- 
cycle of the composite cathode and is fairly consistent across the various RVEs. During the 
second charging half-cycle, no additional damage occurs in any of the simulations and the 
electrode behavior stabilizes. 

In addition, Fig. 3.7(c) shows the evolution in Voltage as a function of the average 
Denem. Plots of this nature are particularly useful in assessing the manner in which adjust- 
ing the voltage caps in which we cycle the electrode can have a significant effect on the 
accumulation of damage. We note again that all RVEs have similar and consistent behav- 
iors. Small differences in behavior of the three RVEs are related to the specific design of 
the microstructure, where for example, interactions between clusters of particles can lead 
to a slightly earlier initiation of damage. However, these differences are not significant 
as compared to changing the volume fraction of active material, as will be discussed in 
Sect. 3.3.2. 

Fig. 3.8 shows the simulation domains at the end of the first half-cycle for the three 
RVEs, where in the active particles we show contours of the normalized concentration, ¢, 
and in the SSE matrix we show contours of the maximum in-plane principal stress, o1. 
Note that deformation is scaled by 10x to best visualize interfacial degradation. 

We observe that damage is consistently distributed across all particles, irrespective of 
particle size, and that we sustain similar levels of stress in the matrix across all simulation 


domains. Finally, these results also hold for microstructure compositions with other active 
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Figure 3.8: Simulation results for composite cathode RVEs with constant active material 
volume fraction dam = 30% and varying particle size distribution. Results are shown at the 
end of the first half-cycle of each simulation domain. Reproduced with permission from 
[17]. 

material volume fractions. We refer the reader to Appendix A.4 for simulation results with 
an active material infill, day = 20%. 

In this Section, we have demonstrated how starting with a fixed composition of ac- 
tive material volume fraction, variations in particle size and particle size distribution do 
not play a significant role in the build up of stress, subsequent damage, and ultimately on 
electrochemical performance. From a design standpoint, this finding has important impli- 
cations as it clearly points to the fact that reducing the size of active particles is not an 
effective mechanism for reducing interfacial stresses and damage. We next discuss the role 


of varying the active material volume fraction, dam, which as expected from the discussion 


surrounding the results shown in Fig. 3.5, can have a significant effect. 


3.3.2 Role of active material volume fraction on electrochemical performance and mechanical 


degradation 


We consider next the role of varying the active material volume fraction, day on elec- 
trochemical performance of SSB composite cathodes. Increasing the amount of active 
material in composite cathodes can be viewed as beneficial for it improves the overall vol- 


umetric capacity of the electrode. However, as demonstrated here, increasing the active 
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material infill can also lead to enhanced stress build up at the active particle/SSE interface, 
which can compromise mechanical integrity. This trade off poses the important question of 
whether there exists a critical active material volume fraction which maximizes volumetric 
capacity in the presence of mechanical degradation at the active particles interface. 

We consider four representative microstructures with active material volume fraction 
in the range dam € [20%, 30%, 40%, 50%] as shown in Fig. 3.9. Particle size distribution 
is kept identical for all simulation domains and is prescribed a uniform distribution with 
particle size lower and upper bounds of 2 and 6 um respectively. As in Sect. 3.3.1 above, 
the aspect ratio of active particles is also uniformly distributed with lower and upper bounds 
set at 0.9 and 1.5 respectively. This size distribution is chosen to match the experimental 
SEM images in the work of Zhang et al. [16]. As before, periodic boundary conditions 
are prescribed to the electrode domain and we cycle the composite cathode over three half- 
cycles between 3.8 - 4.6 V voltage caps at a constant C-Rate of 0.5C. Damage is modeled 


using the same material parameters described in Sect. 3.3.1 above. 
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Figure 3.9: Simulated RVEs with a constant particle size distribution and varying active 
material volume fraction dam. Reproduced with permission from [17]. 


In order to account for increase in capacity with added active material, here, we report 
capacity results normalized by the total mass of electrode. That is, capacity is computed 


as, 
cAVam ` CRmax * F 
Capacity = Irs ee (3.3.1) 


MRVE 


where mryg is the total mass of the composite cathode, Cr max is the maximum molar con- 
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centration of the active material given in Table 3.1, F is the Faraday constant, and c is the 
normalized concentration integrated over the volume of active material. The mass of the 
electrode RVE is computed as Mpve = Vam ` Pam + VssE : Psse, With pam = 4.79 g/cm? 
for LCO as reported in [139, 160], while we determine the density of the SSE, pssg = 2.1 
g/cm? from the work of Zhang et al. [16], who reported both volume and mass ratios for 
various LCO-LGPS microstructure compositions. 

Fig. 3.10(a) and (b) show the Voltage vs. Capacity curves for the four simulation 


domains with and without damage respectively. First, we note as expected, that simula- 
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Figure 3.10: Simulation results for composite cathode RVEs with constant particle size 
distribution and varying volume fraction dam. (a) and (b) show Voltage vs. Capacity (de- 
fined per unit total mass of the RVE) for simulations with and without interfacial damage 
respectively. (c) Stable capacity during cycling as a function of Ham for simulations with 
and without damage and also shows the stable capacity loss due to damage, which is the 
difference of the two simulation results. (d) Voltage as a function of accumulated chemical 
damage (i.e. average Denem in the entire simulation domain) for simulations with varying 
Pam and interfacial damage. Note the legend in (a) applies equally to (b) and (d). Repro- 
duced with permission from [17]. 


tions without damage, Fig. 3.10(b), also show variation in Voltage vs. Capacity response, 
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since capacity is normalized per total mass of the electrode. The effect of damage is 
clearly visible and we observe a drastic change in electrochemical behavior when com- 
paring Figs. 3.10(a) and (b). In particular, in Fig. 3.10(a), we highlight the stable capacity 
of electrode as defined by the stable electrochemical window achieved following the first 
half-cycle. We can clearly observe the significant decrease in capacity when accounting for 
damage at the interface for a given active material infill. 

The stable electrode capacity as a function of the active material volume fraction, dam 
1s shown in Fig. 3.10(c). Here, we plot the stable electrode capacity for simulations with- 
out mechanical damage (green line) and damage at the particle-electrolyte interface active 
(blue line). In addition, we plot the difference in capacity of the two curves and label 
the result as the loss in capacity incurred by mechanical delamination of the interface (red 
line). We observe here the manner in which active material volume fraction plays a critical 
role in electrochemical performance of the electrode in the presence of interfacial damage. 
Comparing the undamaged and damaged curves in Fig. 3.10(c), we observe that as active 
material infill, dam is varied from 20% to 50%, stable capacity shrinks by 31, 46, 57, and 
64 mAh/g. In fact, from Fig. 3.10(c), one can observe that stable capacity remains nearly 
constant with Ham for simulations with active mechanical damage at the interface. This 
1s attributed to capacity losses due to mechanical damage increasing with dam almost as 
rapidly as stable capacity of the undamaged microstructure increases with added active ma- 
terial. That is, any capacity gains achieved by having a higher active volume fraction, Pam 
are almost entirely negated by increasing losses due to mechanical degradation. 

From a design perspective, it is also useful to look at the Voltage vs. Accumulated 
chemical damage behavior of the electrode, as shown in Fig. 3.10(d). Two important trends 
can be observed here. First, damage initiates earlier (i.e. at higher voltages) for microstruc- 
tures with higher dam. Second, the rate of damage evolution (i.e. the slope of the curve at 
intermediate accumulated chemical damage values) is also higher for microstructures with 


higher dam. From a design standpoint, the data shown in Fig. 3.10(d) is also useful in guid- 
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ing one to determine where a voltage cutoff could be placed to reduce mechanical damage 
at interface. For example, at a voltage cutoff of 4.1V, the microstructure with dam = 20% 
will have only incurred approximately 60% accumulated damage as compared to 90% for 
microstructures with dam = 50%. 

The detailed evolution in stress, concentration, and damage across the electrode do- 
main is shown in Fig. 3.11, where we show results at a voltage of (a) 4.5 V (top row), (b) 
4.2 V (middle row), and (c) 3.9 V (bottom row). The corresponding voltages are marked in 


Fig. 3.10(d). In all results, we show contours of the maximum in-plane principal stress, o1 


dam = 20% dam = 30% gam = 40% gam = 50% 
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Figure 3.11: Simulation results for composite cathode RVEs with constant particle size 
distribution and varying volume fraction dam. Results are shown at (a) Voltage of 4.5 V 
(top row), (b) Voltage of 4.2 V (middle row), and (c) Voltage of 3.9 V (bottom row). In all 
results, we show contours of the maximum principal stress, 71, over the SSE matrix, and 
contours of normalized concentration ¢ over the active particles. Note that the colorbar for 
gı is constant over a given Voltage, and the colorbar for ¢ is constant for all simulations 
shown. Deformations are scaled by 10x. Reproduced with permission from [17]. 
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over the SSE matrix, and contours of normalized concentration, ¢ over the active particles. 
In Fig. 3.11(a), at 4.5 V, no damage has yet initiated over the domain (see Fig. 3.10(d)) 
and we can clearly observe the manner in which stresses increase with increase in dam. 
We can also clearly observe how g hotspots arise between nearby particles due to me- 
chanical particle-particle interactions. This trend persists at all voltages. In Fig. 3.11(b), 
at 4.2 V, damage has initiated across all microstructures, with significantly more damage 
occurring across microstructures of higher active material infill, dam. A similar response 
is also observed at 3.9V in Fig. 3.11(c). 

Another important observation is with respect to distribution of normalized concen- 
tration č. Focusing on the results in Fig. 3.11(c) at 3.9V, we can clearly observe that 
concentration of Li is more uniformly distributed over all particles in microstructures with 
small active material volume fraction, Pam. In contrast, simulations with high dam show a 
broad range of concentrations over the various particles. This phenomena arises from the 
presence of mechanical interfacial damage. For simulations with high @am, multiple par- 
ticles have entirely decohered from the surrounding SSE matrix and thus have chemically 
isolated. These particles can no longer store Li-ions and contribute to the overall loss of 
capacity due to mechanical damage experienced by the microstructures with higher active 
material volume fraction. 

In summary, we formulated and numerically implemented a theoretical framework — 
based on use of a chemo-mechanical surface element — for modeling the stress-coupled 
chemo-mechanical interactions in complex electrode designs under galvanostatic charging. 
In particular, the theoretical framework allows one to capture both chemical and mechan- 
ical interactions between active particles in battery electrodes. Under galvanostatic charg- 
ing conditions, the theoretical framework captures chemical interactions between particles 
whereby current is distributed non-uniformly between active particles based on their lo- 
cal stress-coupled interface reaction kinetics. Mechanically, the proposed surface elements 


allow one to capture the role of mechanical interfacial damage on electrochemical perfor- 
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mance by coupling interfacial reaction kinetics (1.e. current) to mechanical damage. This 


allows one to model how loss of mechanical integrity at the interface limits reaction kinet- 


ics and overall electrochemical performance of the electrode. We applied the theoretical 


framework to investigate the chemo-mechanical behavior — including interfacial degrada- 


tion — of a LCO-LGPS composite electrode under galvanostatic charging, with a focus 


on understanding how microstructural features such as local material properties, particle 


size distribution and packing density impact the overall electrochemical behavior of the 


composite electrode. 


The following important findings were made with respect to this electrode: 


We demonstrated that variations in SSE stiffness and active material contraction dur- 
ing lithiation can have a large impact on electrochemical performance of SSB cath- 
odes in the presence of interfacial delamination. The theoretical framework devel- 
oped thus provides an important tool for future exploration of design of composite 
electrodes with different active material-SSE pairs in order to enhance mechanical 


integrity and concurrently electrochemical performance. 


With respect to microstructural composition, we first demonstrated that particle size 
does not play a critical role in stress build up at the interface and subsequent mechan- 
ical degradation. This was demonstrated both through simple single-particle RVEs 


and through complex multi-particle RVEs with different particle size distributions. 


The active material volume fraction was demonstrated to be a critical factor in elec- 
trochemical performance in the presence of mechanical degradation. This microstruc- 
tural descriptor largely dictates the build up of interfacial stresses at the active parti- 


cles, and thus controls the onset and evolution of mechanical interfacial degradation. 


We demonstrated that while increasing the active material volume fraction does nat- 
urally increase the capacity of the LCO-LGPS cathode under consideration, these 


gains can be quickly undermined by the enhanced mechanical interfacial degradation 
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experienced. The theoretical framework developed provides a quantitative analysis 
of the manner in which capacity of a given RVE varies with increase in active mate- 
rial in the presence of mechanical interfacial degradation, and can serve as a useful 


tool for design of complex novel SSB cathode architectures for future applications. 


In chapter 4 to follow, we computationally address another critical issue hindering commer- 
cialization of all-solid-state batteries, namely the phenomenon of Li-metal filament growth 


across solid-state electrolytes with fracture of the latter, causing the battery to short-circuit. 
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CHAPTER 4 
A CONTINUUM ELECTRO-CHEMO-MECHANICAL GRADIENT THEORY 
COUPLED WITH DAMAGE: APPLICATION TO LI-METAL FILAMENT 
GROWTH IN ALL-SOLID-STATE BATTERIES 


We formulate a thermodynamically-consistent, electro-chemo-mechanical gradient theory 
which couples electrochemical reactions with mechanical deformation and damage in solids. 
The framework models both species transport across the solid host due to diffusion/migration 
mechanisms and concurrent electrochemical reaction at crack surfaces within the solid host 
where ionic species are reduced to form a new compound. The theory is fully-coupled 
in nature with electrodeposition impacting mechanical deformation, stress generation and 
subsequent fracture of the solid host. Conversely, electrodeposition kinetics are affected by 
mechanical stresses through a thermodynamically-consistent, physically motivated driv- 
ing force that distinguishes the role of chemical, electrical and mechanical contributions. 
Critically, the framework captures the interplay between growth-induced fracture of the 
solid host and electrodeposition of a new material inside cracks by tracking the damage 
and extent of electrodeposition using separate phase-field variables. Here, the fundamen- 
tal distinction is made such that material deposition across the solid host can only occur 
upon preliminary fracture of the solid host to enable for creation of vacant sites which 
can accommodate the newly deposited compound. An attractive feature of the theory is 
1ts ability to model nucleation, propagation and branching of cracks across the solid host 
in arbitrary orientations driven by a thermodynamically-consistent fracture driving force. 
Critically, the gradient damage formulation invokes a tension-compression asymmetry to 
preserve resistance in compression of fractured regions within the solid host. 

While the framework is general in nature, we specialize it towards a critical problem 


of relevance to commercialization of next-generation all-solid-state batteries, namely the 
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phenomenon of metal filament growth across a solid-state electrolyte, eventually leading to 
short-circuit. We specialize on a Li-metal - Li,LazZr2012 (LLZO) system and demonstrate 
the capacity of the framework to capture both integranular and transgranular crack and 
Li-filament growth mechanisms, both of which have been experimentally observed. De- 
velopment of a microstructure-resolved description of the process of metal filament growth 
through the SSE microstructure is currently lacking and accordingly addressed in this the- 
sis. In modeling this system, we elucidate the manner in which mechanics and fracture 
of the SSE impact electrodeposition kinetics and Li-filament growth. From a manufactur- 
ing standpoint, we additionally elucidate the role of mechanical confinement in solid-state 
batteries on the rate of crack propagation versus the rate of Li-metal deposition across 
the LLZO electrolyte. Under specific mechanical boundary conditions, we demonstrate 
the capacity of the framework to qualitatively reproduce the experimentally observed phe- 
nomenon of crack fronts propagating ahead of Li-metal filaments, as partially-filled cracks 
reach the far electrode end in advance of Li-deposits. Beyond this application, the frame- 
work should serve useful in a number of engineering problems of relevance in which elec- 
trochemical reactions take place within a damage zone, leading to creation of new material 


at these locations. The work in this chapter was adopted with permission from: 


e Bistri, Donald, and Claudio V. Di Leo. ”A continuum electro-chemo-mechanical 
gradient theory coupled with damage: Application to Li-metal filament growth in 


all-solid-state batteries.” Journal of the Mechanics and Physics of Solids 174 (2023). 


4.1 Introduction 


Solid-state batteries (SSBs) present a promising technology for next-generation energy 
storage systems. In recent years, commercialization of SSBs has attracted significant re- 
search attention owing to their promise of superior performance compared to conventional 
Lithium-ion batteries (LIBs). Among their many advantages, SSB architectures can en- 


able for increased current density (~ 3860 mAh/g), improved safety, and a wider electro- 


71 


chemical window (0-5V) - in turn enabling for coupling of Li-metal anodes with high 
voltage cathode materials [3, 4, 5, 7, 161]. From a safety standpoint, use of inorganic 
solid-state electrolytes (SSEs) immediately alleviates potential hazards associated with ig- 
nition of liquid electrolytes commonly used in LIBs. In addition, the solid-state nature 
of the electrolyte material has the potential to alleviate failure mechanisms associated with 
dendrite growth. However, successful design and operation of SSBs is hampered by several 
chemo-mechanical challenges across its constituents [9, 17, 18, 19, 153, 162], the most crit- 
ical one associated with metal filament growth across the SSE microstructure, eventually 
leading to short-circuit. 

Unlike the process of dendrite growth in liquid electrolyte systems, growth of metal 
filaments across the SSE in a solid-state architecture depends on fundamentally different 
mechanics. Owing to the solid-state nature of the electrolyte, the SSE microstructure is 
hypothesized to play a fundamental role. Li-metal filaments can initiate at imperfections 
or heterogeneities of the metal/SSE interface such as voids, cracks, and grain boundaries. 
Once initiated, these filaments propagate through the SSE microstructure by fracturing 
the latter and creating fresh sites for continuing Li-metal deposition and filament growth. 
In this view, mechanical fracture and Li-metal filament growth are intrinsically coupled. 
Eventually, metal filaments may pierce through the entire SSE, reaching the cathode and 
leading to short-circuit of the battery [20, 58, 59, 60]. It is thus critical to understand from 
both an experimental and modeling perspective the interplay between electrodeposition at 
the Li-metal interface and the SSE microstructure. In particular, the manner in which de- 
fects, rate of Li deposition, local stress state, and mechanical behavior (including fracture) 
of both Li-metal and SSE couple to govern the onset and evolution of Li-metal filaments in 
SSB architectures [22, 61, 62, 63, 64, 163]. 

Filament growth has been experimentally observed to occur above a critical current 
density in both compliant [20, 66, 67, 68] and stiff oxide-based [58, 59, 60, 69, 70, 71] 


electrolytes. These experimental observations counter previous modeling predictions [72] 
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that claimed filament growth could be suppressed - and replaced with stable Li-metal depo- 
sition - provided the shear modulus of the SSE is at least twice as large as that of Li-metal. 
Advanced characterization techniques including tomographic imaging and spectroscopy 
have revealed propagation of Li-metal filaments across the SSE microstructure to occur 
either 1) transgranularly, in the form of a dominant Li-filament, typical for single-crystal 
electrolytes [20, 60, 69, 73] or 11) intergranularly, through the grain boundaries, typical for 
polycrystalline SSEs [58, 59, 66, 74, 75]. Recent studies have additionally hypothesized 
the potential for isolated Li-metal deposits to form within the bulk SSE microstructure due 
to the favorable SSE electronic properties and presence of trapped electrons. Upon reac- 
tion with Li-ions, these trapped electrons may produce isolated Li-metal deposits, seeding 
subsequent growth of Li-metal filaments across the SSE microstructure [76, 77]. 

Aside from experimental observations, several modeling frameworks have been pro- 
posed in the literature to address different aspects of the nucleation and propagation of 
Li-metal filaments in SSEs. Related to nucleation of Li-metal filaments, linear perturbation 
analysis has been extensively employed in several models [78, 79, 80, 81, 82, 83]. Start- 
ing with a sinusoidal interface perturbation, these models have investigated the stability of 
electrodeposition at the Li-metal/SSE interface with variation in perturbation amplitude, 
applied current density and state of pre-stress. Consistent with experimental observations, 
growth of Li-metal filaments (e.g. unstable electrodeposition) was shown to occur irrespec- 
tive of the SSE stiffness provided the applied current density and perturbation wavelength 
are sufficiently high. To study the potential for propagation of existing Li-metal filaments, 
a second class of models, treating metal protrusions as pressurized cracks under a linear 
elastic fracture mechanics framework was later developed [60, 164, 165, 166, 167, 168]. 
These works investigated the role that defect size, surface resistance and tip overpoten- 
tial play on the propensity for Li-metal filaments to fracture the SSE and grow across the 
solid conductor. While clearly important, the two classes of models discussed above are 


limited to either modeling the nucleation of Li-filaments or the conditions under which an 
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existing filament might propagate. They can not capture the dynamic evolution of Li-metal 
filaments, which is the focus of this work. 

Recently, mature continuum chemo-mechanical frameworks have been developed to 
model the complex electro-chemo-mechanical processes that govern non-uniform elec- 
trodeposition at a Li-metal/SSE interface in the presence of imperfections. Narayan and 
Anand [84] propose a thermodynamically-consistent continuum theory to model non-uniform 
plating/stripping kinetics at the Li-anode coupled with elastic-viscoplastic deformation of 
Li-metal. Given the complexity of modeling freshly-deposited Li-layers from a finite el- 
ement standpoint, the authors propose an analogous mechanical problem of swelling and 
de-swelling of a thin “interphase layer” at the Li-metal/SSE interface to reproduce the 
electro-chemical process of plating/stripping. The framework provides important insights 
on the role that both geometric imperfections and chemical impurities at the metal/SSE in- 
terface play on the evolution of Li-metal filaments and concomitant deformation and frac- 
ture of the associated metal electrode-solid electrolyte materials. However, the framework 
is limited in its capacity to model the ongoing electrochemical reaction kinetics driving the 
growth of Li-metal filaments, and can not subsequently predict the morphology and evo- 
lution of the resulting Li-filaments. The work of Shishvan et al. [85, 86] presents one of 
the first attempts to model the evolution of Li-metal filaments across the SSE beyond the 
idealization of filaments as pressurized cracks. In [85], the authors propose an alternative 
mechanism for initiation of metal filaments, treating metal protrusions as thick-edge dislo- 
cations and demonstrate the capacity of the framework in [86] to predict the experimentally 
reported critical current density and growth rate of Li-deposits across SSEs. 

Phase-field models have also been extensively developed to capture Li-metal filament 
growth in both liquid and solid-state battery architectures. Purely electro-chemical non- 
linear phase-field models have been developed and shown to successfully predict the evo- 
lution and morphology of Li dendrites in liquid-electrolyte systems under varying applied 


electric potentials [169, 170, 171, 172, 173]. Recently, phase-field based models have 
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been extended to solid-state architectures - including both inorganic and polymer-based 
electrolytes - to incorporate the role of mechanics and deformation of Li-metal due to the 
confining SSE on electrodeposition kinetics [77, 174, 175, 176]. These works present the 
first attempts at capturing the evolution and morphology of Li-metal filaments across a 
solid conductor, while simultaneously accounting for the presence of microstructural inho- 
mogeneities in the solid-electrolyte [174]. However, these models do not capture the fun- 
damental coupling between fracture of the SSE and electrodeposition. Li-metal filaments 
growing in solid electrolytes, unlike dendrites in their liquid counterparts, are constrained 
by a solid continuum and must overcome the associated mechanical resistance by fractur- 
ing the solid host. Only upon fracture of the solid conductor, and subsequent creation of 
an empty deposition site, may Li-metal filaments proceed to grow by plating the newly 
fractured space. Recent work by Fincher et al. [65] experimentally demonstrates this phe- 
nomenon. There, the authors apply an external, experimentally controlled stress field on 
a solid-state electrolyte through which a Li-filament is growing. By controlling the stress 
field, Fincher and co-workers change the preferred fracture direction of the underlying SSE 
and demonstrate that this in turn changes the direction of Li-filament growth correspond- 
ingly. This work provides direct experimental evidence for the coupling between SSE 
fracture and Li-metal filament growth. 

As summarized above, there is a need for a fully-coupled, electro-chemo-mechanical 
framework which captures: i) the nucleation and evolution of Li-metal filaments, ii) the 
concurrent fracture of the solid host, and iii) the coupling between mechanical deforma- 
tion, stress and Li-metal electrodeposition kinetics. Of particular importance in such a 
framework is also the ability to model the role that heterogeneities in the SSE microstruc- 
ture (i.e. pores, cracks, grain boundaries) play on the coupled fracture behavior of the 
SSE and concurrent Li-metal filament growth due to electrodeposition. The purpose of this 
work is to present such a framework. 


We develop here a thermodynamically consistent, phase-field electro-chemo-mechanical 
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framework for modeling concurrent diffusion, reaction, deformation and damage in solids. 


While the framework is general in nature and can be applied across a number of engineer- 


ing problems (i.e. electro-chemically active polymers [57, 177, 178, 179], oxidation and 


corrosion [45, 53, 180, 181, 182], electro-refining [183]), we specialize it towards model- 


ing the onset and evolution of Li-metal filaments in all-solid-state batteries. In particular, 


our framework includes the following unique features distinguishing it from previous work 


in the literature: 


The theoretical framework captures concurrent charged species diffusion/migration 
and reaction, coupled to mechanical deformation and damage of the solid host. Cur- 
rently, continuum electro-chemo-mechanics theories for energy storage materials in 
the literature have been largely developed for diffusion of a conserved species across 
the host (i.e. no chemical reactions) [42, 108, 111, 112, 113], with some recent efforts 
capturing diffusion-reaction in solids [15, 46, 47]. While rigorous, current diffusion- 
reaction-deformation frameworks do not consider transport of charged species across 
the host in the presence of electric field, and thus are limited to modeling chemical 
rather than electro-chemical reactions. Continuum mechanics formulations for trans- 
port of charged species have been developed [87, 177, 178], but in turn do not con- 
sider electrochemical reactions and apply only to a conserved diffusing species. This 
work provides a unified thermodynamically consistent electro-chemo-mechanical 
framework, which concurrently accounts for diffusion and reactions within the host 


material with the potential for reaction-induced fracture of the latter. 


The theory captures the coupling between fracture of the solid host and electrodepo- 
sition of a new material inside cracks by modeling the evolution of cracks and ma- 
terial deposition using distinct phase-field variables. While numerical frameworks 
coupling reaction with fracture of the host material have been proposed in the lit- 
erature (i.e. stress corrosion cracking [52]), the theory importantly presents the first 


framework coupling a phase-field reaction model for Li-metal filament growth with a 
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phase-field damage model for electrodeposition-induced fracture of the solid conduc- 
tor for energy-storage applications. Critically, this enables one to model the manner 
in which mechanical fracture of the SSE is intrinsically coupled to electrodeposi- 
tion and growth of Li-metal filaments. We demonstrate how this coupling can arise 
both through a thermodynamically consistent reaction driving force, and/or through 


a damage dependent formulation of the reaction kinetics. 


The theoretical framework, and accompanying numerical implementation, allows 
one to model the role of microstructure heterogeneities (i.e. voids, cracks, grain 
boundaries) with varying chemo-mechanical properties on concurrent fracture of the 
SSE and electrochemical growth of Li-filaments. As such, the model may predict the 
manner in which various microstructural features enhance or suppress Li-filament 
growth and subsequently guide the development of microstructures for minimizing 


filament growth-induced degradation. 


To demonstrate the relevance and use of the proposed theoretical framework, we spe- 
cialize it towards modeling the growth of Li-metal filaments in an inorganic Lithium Lan- 
thanum Zirconium Oxide electrolyte, specifically Li7La3Zr2O012 (LLZO). LLZO is a promis- 
ing SSE candidate due to its high ionic conductivity, high mechanical stiffness and good 
stability against Li-metal. We explore the experimentally observed manner in which Li- 
metal filaments may grow across the LLZO electrolyte either 1) transgranularly, as a single 
dominant Li-metal filament propagating across the conductor [60], and/or 11) intergranu- 
larly, across the mechanically weaker grain boundaries [58, 59, 66]. In particular, owing 
to the coupling between mechanical stresses, fracture and electrodeposition kinetics devel- 
oped in this framework, we investigate the manner in which mechanically weaker grain 
boundaries with lower fracture energy dictate the transgranular versus intergranular na- 
ture of the resulting Li-metal filaments. Finally, using our framework, we provide insights 
on the manner in which mechanical confinement and stresses dictate the advancement of 


cracks versus Li-metal deposits across the solid electrolyte. In particular, whether the crack 
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front propagates significantly ahead of the Li-metal filaments as has been experimentally 
observed [20, 66]. 

The framework is organized as follows. In Sect. 4.2, we introduce mass balance, charge 
balance and electrostatics equations and define the physically motivated phase-field param- 
eters, €, governing the extent of electrodeposition, and d, governing the extent of damage 
in the solid host. Kinematics are developed in Sect. 4.3. The remaining governing laws are 
developed in Sect. 4.4 using the principle of virtual power and the first and second laws 
of thermodynamics. The constitutive theory is presented in Sect. 4.5 and summarized in 
its general form in Sect. 4.6. In Sect. 4.7, we present a specialization of the theory to- 
wards modeling the growth of Li-metal filaments in a LLZO solid-state electrolyte, and 
summarize the governing partial differential equations and initial/boundary conditions in 
Sect. 4.8. Numerical simulations are presented in Sect. 4.9. First, in Sect. 4.9.1, we model 
the growth of a single Li-metal filament across a LLZO electrolyte. Through this model, 
we demonstrate the salient features of our theoretical framework and explore the role that 
mechanical confinement and stress play on the evolution of cracks and Li-metal filament 
growth across the solid electrolyte. In Sect. 4.9.2, we illustrate the manner in which dif- 
ferent microstructural features, in particular grain boundaries, may be incorporated into 
the modeling framework. Within this model of a LLZO electrolyte with resolved grain 
boundaries, we elucidate the manner in which mechanical stresses, damage and electrode- 
position are coupled and how different coupling mechanisms result in different fracture 
and Li-filament morphologies. Finally, we demonstrate the manner in which this frame- 
work directly couples mechanical properties, such as grain boundary fracture energy, with 


the resulting Li-filaments growth pattern. We close with concluding remarks in Sect. 4.10. 


4.2 Phase-field formulation and balance equations 


This section details the conservation laws for electrodeposition phenomena within a solid 


ion-conducting host. The problem is shown schematically in Fig. 4.1(a), where we illus- 
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trate a solid conductor (SSE) adjacent to an electrode composed of a metallic compound 
“M”. The solid conductor/metal electrode interface may contain a number of defects (i.e. 
pores, voids, cracks), which serve as vacant sites for the metallic compound to fill and 
subsequently deposit inside the solid host (c.f. [60, 164, 165]). Consider now a simple 
and general electrodeposition reaction. Cations, M”* conduct across the solid host towards 
the conductor/metal interface, where they react with electrons, e” to deposit M-atoms as 
follows 


M"t + ne” > M. (4.2.1) 


The electrodeposition of M-atoms (i.e. metallic compound) is treated in a phase-field sense 
and taken to occur over a diffuse boundary, with the extent of reaction tracked via a nor- 


malized phase-field parameter 


E =E/E"" € [0,1]. (4.2.2) 


Here, € denotes the moles of electrodeposited species per unit reference volume, with €™** 
the maximum amount of the metallic compound that can be electrodeposited at the reaction 
site. Physically, a state of € = 0 denotes the absence of metallic compound or “void”, while 
a state of £ = 1 denotes the fully deposited metallic compound. Naturally then, intermediate 
values of 0 < € < 1 represent the reaction zone. 

In this physical interpretation of the phase field variable, €, we restrict the formulation 
such that electrodeposition inside the solid host may occur only within voids (i.e. cracks, 
pores), which either form through mechanical loading or exist inherently in the microstruc- 
ture. This interpretation conforms with experimental observations that for metal deposits 
to grow through a solid host, they must first overcome mechanical resistance by fracturing 
the solid to create the necessary vacant space to accommodate plating (c.f. Fincher et al. 


[65] and Ren et al. [184]). 


We now introduce a phase-field damage variable, d(X,t) to describe any state of im- 
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perfections across the solid host, either pre-existing (i.e. pores, voids) or mechanically 
induced (crack formation) 


d(X,t) e [0, 1]. 


Here, a state of d = 0 represents the intact material, while a state of d = 1 represents the 
fully-fractured solid host. Intermediate values, 0 < d < 1 naturally denote then a partially 
fractured material. Additionally, we assume microstructural changes leading to fracture to 


be irreversible and require damage to grow monotonically such that 


d(X,t) > 0. (4.2.3) 


With a focus on the crack/electrodeposition interplay, we restrict ourselves such that 


electrodeposition may only occur within damaged regions of the SSE. That is, 


E>0 onlyif d>0. (4.2.4) 


This is shown schematically in Fig. 4.1(b-d). Fig. 4.1(b) illustrates a microstructural defect 
in the form of a void or crack (region of d = 1) at the metal/SSE interface, which provides 
the necessary vacant space for metal to deposit. With continuous deposition, the metallic 
compound progressively covers the entire area of the crack (d = 1, € = 1), as illustrated 
in Fig. 4.1(c). This in turn induces a build-up in stresses acting on the crack’s surface, 
large enough to eventually overcome the fracture toughness of the solid conductor and 
incur further damage on the solid host, Fig. 4.1(d). The newly formed damaged region 
(d = 1,€ = 0) then provides the additional vacant site for the metal to deposit further into 
the solid electrolyte. This synergistic mechanism sustains the growth of Li-metal filament 
across the solid host. Eventually, the metallic protrusions fracture and transverse the entire 


span of the solid conductor, at which point the battery short-circuits. 
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Figure 4.1: (a) Schematic of a metal electrode-solid conductor assembly. (b) Illustrates 
a pre-existing defect (€ = 0,d = 1) at the solid host-metal electrode interface prior to 
metal deposition inside the crack. (c) Illustrates deposition of metal inside the imperfection 
(€ = 1,d = 1), leading in turn to a build up in normal stresses, cp, at the crack surfaces. 
(d) Build-up in normal stresses on the crack surface eventually causes the solid host to 
damage further ahead of the crack tip, creating a new damage zone (€ = 0,d = 1), which 
can accommodate subsequent deposition of metal inside the solid host. Reproduced with 
permission from [19]. 


4.2.1 Mass Balance 


Considering the electrodeposition reaction (4.2.1), conservation of mass for the mobile 
ionic species transporting across the solid host may then be written as a diffusion-reaction 
equation of the form 


è= —Divj, — Ê. (4.2.5) 


Here, c(X,t) denotes the number of moles of the mobile ionic species per unit reference 
volume, while £(X, t) denotes the number of moles of the electrodeposited species per unit 
reference volume with € the electrodeposition rate. Consistent with the discussion in Sect. 
4.2, when the extent of reaction E = 1, the reaction has led to consumption of €™** moles 


of the diffusing species. Additionally, jr (X, t) represents the referential flux of the mobile 
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ionic species transporting across the solid host, which must be constitutively prescribed. 
Note here that the reaction-diffusion equation (4.2.5) encapsulates both transport of ionic 
species across the solid host via the Div jr term and its subsequent consumption leading to 


formation of the newly electrodeposited metallic compound via the € term. 


4.2.2 Charge Balance 


The net charge per unit reference volume of the solid conducting host, q, is given by (c.f. 


Narayan and Anand [177] and Li and Monroe. [185]) 


de = F(c— co), (4.2.6) 


where co represents the initial (i.e. bulk) concentration per unit reference volume of the 
mobile ionic species transporting across the solid host, and F is the Faraday constant. Note 
here that in writing Eqn. (4.2.6), we restrict our attention to monovalent cationic species 
and additionally assume negatively charged species in the solid host to be immobile. This 
restricts the mass flux across the solid host to be carried only by the mobile cations. As a 
result, any local excess charge density across the solid host can be induced only by a change 
in concentration of the mobile cationic species from the equilibrium bulk concentration as 
they transport across the solid conductor. 


Defining a referential current density, 


i Pjk (4.2.7) 


taking the time derivative of (4.2.6), and making use of mass balance (4.2.5), one may write 
for later use 


dk = —Divi, — FE. (4.2.8) 
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4.2.3 Electrostatics 


Under electrostatic conditions, for a polarizable material we require that the two governing 
Maxwell equations are satisfied. The first equation, Faraday’s law, requires the referential 


electric field, e,(X, t) to obey the following relationship 


curle, = 0. (4.2.9) 


One can automatically satisfy Eqn. (4.2.9) by representing the electric field, e,(X, t) as the 


gradient of an electrostatic potential, p(X, t) such that 


e,(X,t) = —V4(X, t). (4.2.10) 


The second equation of electrostatics is given by Gauss’s law, and relates the electric dis- 


placement, d,(X, t) to the net charge per unit reference volume, qr as follows 


Div dy = qr. (4.2.11) 


In conjunction, eqs. (4.2.9) and (4.2.11) constitute the governing equations for electrostat- 


ics in the solid host. 


4.3 Kinematics 


We introduce here the kinematic formulation for description of a deformation resulting 
from coupled species transport, distortion due to electrodeposition of a new material inside 
the solid host, and mechanics. Consider a macroscopically homogeneous body B with the 
region of space it occupies in a fixed reference configuration and let X denote an arbitrary 


material point of 5. The motion of B is then a smooth one-to-one mapping x = x(X, t) 
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with deformation gradient, velocity, and velocity gradient given by*. 


F= VX, v=x, L = grad v = FF™!. (4.3.1) 


The theory builds on a multiplicative decomposition of the deformation gradient, 


F — Frmechanical Frchemical — PF FR: (4.3.2) 


Here 


(i) F°(X) represents the concurrent local distortion of the material neighborhood of 
X due to chemical phenomena including: i) species transport across the solid host 
electrolyte, and ii) distortion due to electrodeposition of a new material inside the 


solid host at the reaction sites. 


(ii) F™ represents the distortion due to macroscopic stresses and may include local irre- 
versible plastic deformation caused by inelastic mechanisms such as dislocation mo- 
tion, and the subsequent elastic stretching and rotation of the inelastically deformed 


material neighborhood. 


Remark 1: In modeling local distortions due to electrodeposition through the tensor 
F°, the chemical deformation tensor acts as a surrogate for modeling the continuous 
creation of new material and ensuing distortions, without explicitly accounting for the 
temporal evolution of the domain. Similar continuum treatments with a focus on growth 


of electrochemical interphases and growing matter in living systems can be found in the 


works of Narayan and Anand [84], Tantratian et al. [174] and Kuhl [187]. 


“Notation: We use standard notation of modern continuum mechanics (c.f. Gurtin et al. [186]). Partic- 
ularly: V and Div denote the gradient and divergence operators with respect to the material point X in the 
reference configuration; while grad and div operate on the point x = x(X, t) in the deformed body; a super- 
posed dot denotes the material time-derivative. Throughout, we write F°-! = (F*)71, F*7 = (F*)7, etc. 
We also write tr A, sym A, skw A, Ao, and sym,A respectively, for the trace, symmetric, skew, deviatoric, 
and symmetric-deviatoric parts of a tensor A. Finally, the inner product of tensors A and B is denoted by 
A: B, and the magnitude of A by |A| = VA: A 
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We refer to F™ and F° as the mechanical and chemical distortions respectively. The 
volume ratio is given by 


JE detF > 0, (4.3.3) 


and using (4.3.2), 


J=J®J°, with J™“detF™>0, and JJ: detF° > 0, (4.3.4) 


such that F™ and F° are invertible. The right and left polar decomposition of F™ is given 
by 
PS RUS VR (4.3.5) 


where R™ is a rotation, while U™ and V™ are symmetric, positive-definite right and left 
stretch tensors. Consistent with standard notation, the total and the mechanical right Cauchy- 


Green deformation tensors are given by 


C=FE and C™ = F™ F". (4.3.6) 


Using (4.3.1) and (4.3.2), the velocity gradient may be written as 


L = L” LE (4.3.7) 


with 


L” = FF"! L' = FF, (4.3.8) 


We define the mechanical and chemical stretching and spin tensors through 


D” =symL”, W” = skw L”, 
(4.3.9) 
D° = symL*, W° = skw L°, 
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such that L” = D™ + W™ and L° = D° + W°. 
Following Afshar and Di Leo [15], we employ a weighted decomposition of the chem- 


ical velocity gradient of the form 


L° = (1—h(£))L* + h(£)L" (4.3.10) 


with h(€), an interpolation function of the reaction coordinate £, which satisfies h(0) = 0 
and h(1) = 1. This function is specified in Sect. 4.7 and consistently used for all quantities 
interpolated by the phase-field parameter €. 

The weighted decomposition (4.3.10), captures the combined deformation due to species 
transport across the solid host and distortion due to electrodeposition of a new material at 


the reaction sites. Here: 


i) L* captures the deformation of a material point due to transport of ionic species 


across the solid host; 


11) L" captures the deformation due to electrodeposition of a new material inside the 


solid host. 


Scaling by the interpolation function h(€) ensures that chemical distortions at material 
points in the bulk of the solid host, away form the reaction sites (i.e where € = 0) are only 
due to transport of ionic species, while in the fully-reacted state (i.e. where € = 1) due to 
electrodeposition of a new compound. 

While the discussion on kinematics is so far presented in its most general form, we spe- 
cialize our work specifically to modeling single-ion conductors. This class of conductors 
is restricted to a single mobile cationic species transporting across the solid host, while 
anions remain immobilized in the backbone of the solid conductor. In such a material, both 
temporal and spatial concentration gradients are negligible within the bulk of the solid host 
(i.e. away from the reaction sites), where deviations from electroneutrality (i.e. equilibrium 


bulk concentration) are insignificant [188, 189, 190]. Additionally, as detailed in [79, 164, 
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191, 192], molar volume of Li-ions within the solid host electrolyte vanishes based on the 
notion that Li-ions lie within a rigid ceramic skeleton of the electrolyte that does not de- 
form upon removal/addition of a Li atom. Henceforth, no deformations are associated with 


transport of ionic species across this class of solid ionic conductors, and hence LY = 0. 


Remark 2: The choice to present the kinematic formulation in its most generic form is 
intended to showcase the versatility of the framework through the kinematic decompo- 
sition (4.3.10) to model the phenomena of metal filament growth for the more general 
class of binary solid-ionic conductors (i.e. polymer solid electrolytes). For this class of 
solid conductors, chemical deformations may arise due to both a combination of species 


diffusion across the host and distortion due to electrodeposition of a new material at the 


reaction sites (c.f. Ganser et al. [193]). 


Focusing now on single-ion conductors, deformations associated with transport of ionic 


species across the solid conductor can be neglected, and (4.3.10) simplifies to 


L° = (8 L”. (4.3.11) 


Further, we make the pragmatic assumption that chemical deformation is irrotational (i.e. 
W° = 0) so that 


D° = h(£)D', (4.3.12) 


with D", the electrodeposition-induced stretching assumed to depend on the reaction rate, 
€ through 
D' = €N' (4.3.13) 


where N' denotes the direction of electrodeposition-induced deformations. 
Finally, the formulation presented so far makes no assumption regarding the nature of 
F™, which may be further decomposed to model, for example, the elastic-plastic behavior 


of the underlying solid conductor host and/or of the newly electrodeposited metal. We now 
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restrict our formulation such that all mechanical deformations are purely elastic in nature 
and further specialize F™ = F°, with F* the elastic mechanical distortion. That is, both the 
solid conductor host and the newly electrodeposited metal will behave elastically. 

Finally, using (4.3.1), (4.3.2), (4.3.11) and (4.3.13) we may then rewrite (4.3.7) for 
future use as 


(Vx)F-! = FF + A(€E)F°(EN) EO. (4.3.14) 


Remark 3: With regards to the solid ionic conductor, the assumption of purely elastic 
deformations is consistent with their stiff brittle nature. On the other hand, continuous 
deposition of metallic material (e.g. Li-metal) inside defects induces large compressive 
stresses, which can cause the metal to flow plastically. As the metal flows plastically, 
it may extrude from the crack to the electrode, limiting the pressure that can build up 
inside the crack and in turn the filament growth behavior. 

As discussed by Klinsmann et al. [164] and Bucci and Christensen [165], the possi- 
bility for the metal to plastically flow out of the crack mouth is impractical for realistic 
battery designs, suggesting that metal remains confined within the crack. In such a 
case, pressure will build up in the metal-filled defect until further damage occurs or the 
reaction is suppressed by mechanical stresses, causing the metal deposition to cease. 
The nature and relevant importance of Li-metal plastic flow on filament growth through 
solid-state electrolytes remains a topic of active research. In this thesis, we restrict our 
theoretical formulation to modeling Li-metal deformations as purely elastic. Owing to 


the already complex phenomena being modeled in this work, we purposely leave the 


addition of modeling an elastic-plastic metal phase to a future publication. 


4.4 Governing Balance Laws 


In this section, we develop the governing equations for our theoretical framework, includ- 


ing macroscopic and microscopic force balances and thermodynamic laws. For concise- 
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ness, we relegate portions of the detailed development to Appendix B, and present here the 


critical aspects of developing the governing balance laws. 


4.4.1 Principle of virtual power. Balance of forces 


To develop the remaining balance laws in our theoretical framework, we invoke the virtual- 
power approach (cf. Gurtin et al. [186]). This results in a macroforce balance and mi- 
croforce balances for the rate-like kinematical descriptors in our theory. In exploiting the 
principle of virtual power, we consider a list of generalized virtual velocity fields to be 
given by 

V = (9x,0F*, 0€, VÓ£, ôd, Vdd). (4.4.1) 


In light of (4.3.14), the virtual velocities are not independent and must obey the following 
kinematic constraint 


(Vóx)F"? = 6F°F*! + h(€)F°(6EN') Fo (4.4.2) 


For any part P of the reference body B, the internal and external power are formulated in a 


consistent fashion with the manner in which forces expand power on P, 


Weal P, V) = ll 


ta(n) dx dos + [by dx dun + | meda, + | yôd dag 
OP P OP OP 


(4.4.3) 
Win P, V) = [os F° + ESE + G- VOE + wôd + Ç: Vdd) dug 
P 


where we have defined the following macroscopic and microscopic force systems conjugate 


to the kinematical variables, 
a) A stress S° that expends power over the elastic distortion rate F°; 
b) A traction t,(n,) (for each unit vector nz) that expends power over the velocity X; 


c) A scalar microscopic stress E that expends power over the rate E; 
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d) A vector microscopic stress G that expends power over the gradient VE 


e) A scalar microscopic traction 7(n,) that expends power over € on the boundary of 


the part; 
f) A scalar microscopic stress w that expands power over the rate d; 
g) A vector microscopic stress ¢ that expands power over the gradient Vd; and 


h) A scalar microscopic traction y(n,) that expands power over d on the boundary of 


the part. 


Note that consistent with the discussion in Sect. 4.3, no deformations are associated with 
the transport of ionic species in single-ion conductors. As a result, no mechanical power is 
associated with the concentration variable c. 


The principle of virtual power consists of two basic requirements: 


i) Power balance, which requires that 9Wx(P, V) = dWin(P, V) for all generalized 


virtual velocities Y 


ii) Frame-indifference, which requires that )Win(P, V) is invariant under all changes in 


frame. 


We relegate a detailed development of the derivation of macroforce and microforce bal- 
ances for the rate-like kinematic descriptors in our theory to Appendix B. Here, we suc- 
cinctly summarize the macroforce balance for the Piola stress, T} as well as the two corre- 


sponding microforce balances for the stresses { E, G, ¢, w}. These are given by 


Div Tk + bk = 0, with boundary condition t,(ng) = Tang, 


E — J°h(€)M*:N' — DivG = 0, with boundary condition (ng) = G- ng, 


Div¢ —@=0, with boundary condition y(nz) = ¢- Ne. 
(4.4.4) 
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In writing (4.4.4), we have introduced the classical Piola stress Tk = JTF~', with T the 


Cauchy stress. In addition, we define the elastic Mandel stress as 


Me jeep”, (4.4.5) 


The macro and microforce balances, when supplemented with a set of thermodynamically- 
consistent constitutive equations, provide the governing mechanical equations for the the- 
ory. We note here that a discussion on invariance properties of the fields, consequence of 
the invariance of internal power under a change in frame, is omitted here for the sake of 
briefness. We refer the reader to the works of Anand [113] and Afshar and Di Leo [15] for 


detailed derivations in the context of similar diffusion-reaction-deformation frameworks. 


4.4.2 Balance of energy. Entropy imbalance. Free energy imbalance 


Our discussion of thermodynamics involves the following fields 


Ek the internal energy density per unit reference volume, 

Mk the entropy density per unit reference volume, 

dx the heat flux per unit reference area, 

de the external heat supply per unit reference volume, 

Y the absolute temperature (Y > 0), 

H the chemical potential, 

@ the electric potential, 
and follows the discussion of Gurtin et al. [186]. Consider a material region P and assume 
inertial effects (i.e. kinetic energy) to be negligible. Under isothermal conditions, the two 
laws of thermodynamics reduce to a single statement requiring that the rate of increase 


in free energy of any part P is less than or equal to the power expended on P. Letting 


Vr denote the free energy per unit reference volume, this requirement takes the form of a 
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free-energy imbalance, 


Jo. dur < Wext(P) == I TAN : Mz dar z f Fojr : Mz dag = ¿de : Mz dag (4.4.6) 
P OP OP OP 


Here, the second term on the right hand side of (4.4.6) represents the flux of energy carried 
into P by the flux, jr of the ionic species transporting across the solid host. Additionally, 
following Kovetz [194], the power expended on P by external charges is represented by 
an electromagnetic energy flux, which in the electrostatic limit, is given by the last two 
boundary integral terms on the right hand side of (4.4.6). 


Consistent with electrochemical notation, we define 


e def 


u+ Fo (4.4.7) 


to denote the electrochemical potential of the mobile ionic species transporting across the 
solid host, with F the Faraday constant. Making use of (4.4.7), one can alternatively write 


the energy imbalance (4.4.6) as 
[vs dur < Wext(P) — I [jx ` Nrdar — ódz - Ny dar. (4.4.8) 
P OP OP 
Defining the elastic second Piola stress as, 
Te E yepe-lppet = Cele (4.4.9) 


with C° = F°'F* the elastic right Cauchy-Green tensor, the stress power S° : F° in (4.4.3) 
may be written as S°: F° = (1/2) J“T*: C°. Equating the external power with the internal 


power (4.4.3), and applying the divergence theorem to the boundary integral terms, the 
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energy imbalance (4.4.8) may be written as 


1 (ve T an C° — EE — G ; VE — wd — Ç : Vd + u°Div jr + ¢Div dgr 
p (4.4.10) 
Osa EE E Vp) dur <0 
Making use of mass balance (4.2.5), charge balance (4.2.6), Gauss’s law (4.2.11), along 
with the definition of electric field (4.2.10) and electrochemical potential (4.4.7), and using 
the fact that this inequality must hold for all parts P, we obtain the local form of the free 


energy imbalance as 
De— GUT: OC — EEG Vg-wd-Ç: Vd- pt- uE — er: dr +j: Vu? <0. (4.4.11) 


For later use, we define the dissipation D > 0 per unit volume per unit time as 


D = SIT: C+ EE+G-VE+ md +6 Vat pet ¿er deje Vu Ux > 0. (4.4.12) 


Finally, we note here that all quantities in the free energy imbalance (4.4.11) and dissipation 
inequality (4.4.12) are invariant under a change in frame (c.f. Anand [113] and Afshar and 


Di Leo [15]). 


4.5 Constitutive Theory 


We divide this section into energetic and dissipative constitutive equations. 


4.5.1 Energetic Constitutive Equations 


Guided by the free-energy imbalance (4.4.11), we consider a set of constitutive equations 


for the free ener r the stress T°, the chemical potential u, the vector microforce € and 
8y p H 
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the electric field e, of the form 
ba =Úr (A), T"=T (A), w=A(A), C=C(A), and ex=éx(A), (4.5.1) 


where A denotes the list 


A = (C°, c, £, VE, dg, d, Vd). (4.5.2) 


Substituting the constitutive equations (4.5.1) in the free-energy imbalance (4.4.11) yields 


OW 1 cre |. Gre | OW, : On, elé¢ Ove : 
(Se -jrr'):¢ k ne E E ex (3 c) VE 
be _ { Ode [ve Pim Ss 
NE -waat (5% —¢) vl ($e) sh ~ mad iV < 0, 


(4.5.3) 


where we have introduced a decomposition for the scalar damage microstress of the form 


W = Wais + Wen- (4.5.4) 


Motivated by Anand et al. [48], the decomposition (4.5.4) allows for both energetic and 
dissipative effects associated with temporal changes in d, while effects due to the gradient, 
Vd, a measure of the inhomogeneity of damage at the microscale, are taken to be entirely 
energetic. 

From an electrochemical standpoint, in (4.5.3), we make the distinction that processes 
associated with diffusion (governed by ¿) and migration (governed by dg) are energetic, 
while the ones associated with electrodeposition (governed by €) are dissipative. An ex- 
ception to such choice is the power conjugate to Vé, which is taken to be entirely energetic 
(i.e all reaction dissipative processes are already accounted for in the term é). 

As the inequality in (4.5.3) is to hold for all values of [C*, ċ, VË, d, Vd, dp}, their 


“coefficients” must vanish, for otherwise they may be chosen to violate (4.5.3). We are 
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therefore led to the thermodynamic restriction that the free energy determines the stress 
T°, the chemical potential u, the vector microstress G, the scalar microstress Wen, the 


vector microstress Ç and the electric field e, through the “state relations” 


OLEA) 

. ace ” 
_ OYR(A) 
7 dc i 
MENA 
OVE * 
= Otie(A) 
Wen P od ) 
qa Aprl A) 
ovd i 

= Ppl A) 

Od, 


eee hae 


(4.5.5) 


er 


Recalling (4.4.7) and using (4.5.5)2, the electrochemical potential of the ionic species trans- 


porting across the solid host, u° is given by 


pe =) (4.5.6) 


C 


4.5.2 Dissipative Constitutive Equations 


The reduced dissipation inequality is now given by 


-N (ae 


oe Bes r) E+ waisd — jn Vue > 0. (4.5.7) 


We define the electrochemical potential of the electrodeposited species as 


us def OEA) 


a E. (4.5.8) 


Then, consistent with the notion of electrochemical reactions, we define the thermodynamic 


driving force for electrodeposition as a difference in electrochemical potential of the species 
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participating in the reaction through 


def 


FE pwp. (4.5.9) 


Using (4.5.9) in (4.5.7), the dissipation inequality may be rewritten as, 


D = —FE + wad + Vu? - Minos Ve > 0, (4.5.10) 


where guided by the dissipation inequality, we introduce a Fick-type relation for the flux, 


jr Of ionic species across the solid host of the form 


Je = -Mmo V u° (4.5.11) 


where as is standard, the flux, j, is linearly proportional to V u° through the mobility tensor 
M mob- 


Finally, we assume all terms in (4.5.10) to individually satisfy, 


—FE > 0, 
Waisd > 0, (4.5.12) 


Vu : Mmo V u? > 0. 


With respect to dissipation due to electrodeposition (4.5.12), we assume, consistent with 
electro-kinetics theory, that € > 0 if and only if F < 0, and vice versa € < 0 if and only 
if F > 0 (c.f. Bazant [126], Fuller et al. [189]). Given the constraint (4.2.3) for damage 
irreversibility such that d > 0, we satisfy (4.5.12). by requiring wais > 0. Finally, (4.5.12) 


leads to the restriction that the mobility tensor Mmob is positive semi-definite. Altogether, 


these restrictions ensure that the dissipation inequality (4.5.10) is not violated. 
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4.6 Summary of the general constitutive theory 


In this section, we summarize our general phase-field diffusion-reaction-damage theory. 


The theory relates the following fields: 


x= x(X, t), motion; 

F= Vx, J=detF >0, deformation gradient; 

F = EYES, multiplicative decomposition of F; 
F°, J° = det F° > 0; elastic distortion; 

F°, J? = det F° >Q, chemical distortion; 

L = FF! = L° + F*L*F*. velocity gradient; 


L° = FF! = D° + W° with W°=0 chemical velocity gradient; 


D° = h(€)D* chemical stretching; 

D' = €N’, electrodeposition induced stretching; 

F° = R°U® = VR, polar decompositions of F®; 

C° = (F°) F° = USA right Cauchy-Green tensor; 

T= Cauchy stress; 

M° = J'E TES,” elastic Mandel stress; 

T= JTE"; Piola stress; 

T° = JE TE, elastic second Piola stress; 

Ur, free energy density/unit reference volume; 

6; moles of diffusing species/unit reference volume; 
É, moles of deposited species/unit reference volume; 
E = E/E e [0,1], extent of the reaction; 

VE, gradient of the reacted species concentration; 

es electrochemical potential of diffusing species; 
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des referential species flux vector; 


ep, referential electric field; 

Q, electric potential; 

dx referential electric displacement; 
d phase-field damage variable; 

Vd gradient of the damage variable. 


4.6.1 Kinematics and Free Energy 


The deformation gradient is decomposed as 
F =F*F", (4.6.1) 


with F* the elastic distortion, and F° the chemical distortion. Further, the chemical velocity 


gradient is given by L° = F*F“-! = D°, with the chemical stretching D° specialized as 
D° =h(€)D' and D'=ÉN. (4.6.2) 


Here, D" is the electrodeposition induced stretching, with " the reaction rate, N' the direc- 
p 8 


tion of electrodeposition induced deformations, and h(£), an interpolation function, which 


satisfies h(0) = 0 and h(1) = 1. The free energy is given by 
Y =b(A) with A thelist A = (C°, c, £, VE, de, d, Vd). (4.6.3) 


With the direction of electrodeposition induced deformations, N" and the free energy func- 


tion, (A) prescribed, the following quantities may be derived. 
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4.6.2 Stress. Force balance 


The second Piola stress is given by (4.5.5); 


OLLA) 


TOS hoe 
aCe 


(4.6.4) 


with the Cauchy stress given by T = J°-'F°T°F®", the elastic Mandel stress by M° = 
JE" "TE", and the first Piola stress by Tk = JTF ™'. 


The stress is governed by force balance (4.4.4), viz. 


Div T, + br = 0, (4.6.5) 


with Neumann boundary condition (4.4.4);, t.(nz) = TM, where tp is a prescribed trac- 


tion. 


4.6.3 Electrochemical potential. Flux. Mass Balance 


Using (4.5.5) in (4.4.7) and recalling (4.5.11), the electrochemical potential and referential 


flux for the mobile ionic species are given by 


Prl A : 
ve = onl) + Fo and Jr = —M mob V u° (4.6.6) 


with Mob, a positive semi-definite mobility tenor. 
The concentration of ionic species transporting across the solid host, c is governed by 
mass balance (4.2.5), viz. 


è= —Divjk — Ê. (4.6.7) 


with €, the species electrodeposition rate at the reaction site. 
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4.6.4 Reaction driving force. Electrodeposition Kinetics 


Combining (4.5.8), (4.5.9), (4.4.4) and (4.5.5)3 yields the electrodeposition driving force 


as 


F=- pu, with p= GU J*h(E)M*:N' — Div (2) (4.6.8) 


The reaction kinetics are constrained by the dissipation inequality (4.5.12), to obey — Fé > 
0. Note that due to the phase-field (gradient) nature of the extent of reaction £, (4.6.8) 


constitutes a PDE. 


4.6.5 Electrostatics 


The electric potential @ is governed by Gauss’s Law. Combining (4.2.11) with (4.2.6) and 
(4.5.5)¢ with (4.2.10) yields 


Div (dk) = F(c—c), with e, = and e, = —V¢(X,t). (4.6.9) 


4.6.6 Damage 


Combining (4.4.4)3 with (4.5.5)4 5 and recalling the decomposition of the microstress w = 
Wdis + Wen from (4.5.4) yields the governing equation for the phase-field damage variable 


as, 


— ( ab(A)\ _ A) = 
Div ( avd Wais = 0, (4.6.10) 


where the constraint Wais > O is required by the dissipation inequality (4.5.12)». 


4.7 Specialization of the constitutive equations 


The theory developed above and summarized in Sect. 4.6 is quite general. We introduce 


now a Set of specialized constitutive equations to elucidate: i) the process of plating and 


100 


fracture due to electrodeposition of Li-metal in single-ion conductors, specifically a LLZO 
inorganic electrolyte, and ii) the inherent coupling of electric, chemical and mechanical 
effects on electrodeposition kinetics and evolution morphology of Li-metal filaments across 
the solid-state electrolyte. 

The general electrodeposition reaction (4.2.1) is now specialized for Li-metal plating 
as 


Li? +e > Li. (4.7.1) 


4.7.1 Electrodeposition-Induced Deformation 


We begin by specifying the manner in which electrodeposition of Li-metal at the reac- 
tion sites induces mechanical deformations, specifically the direction tensor N" in (4.3.13). 
Consistent with Narayan and Anand [84] and Rejovitzky et. al [44], we consider a generic 


distortion due to electrodeposition such that 


F' = DE (4.7.2) 


Additionally, we assume the electrodeposition induced stretching tensor D" is anisotropic 
and specialize it as 


D'=€S' where SEm,om. (4.7.3) 


Here, the unit vector, m, represents the direction of electrodeposition induced deforma- 
tions, while €" = tr(D") denotes the total volumetric strain rate for reaction induced defor- 
mations. Consistent with experimental observations [195], we assume deformations due to 
electrodeposition of Li-metal to occur preferentially in a direction normal to the reaction 


front (e.g. normal to the conductor/metal interface), and hence specialize m, as 


m, = VE/|VE]l. (4.7.4) 
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Further, we assume volumetric distortions due to electrodeposition to vary linearly with 
the extent of reaction (c.f [113, 177]), such that JT = 1 + QE. Then, using the fact that 


J’ J‘! = tr(D‘), and combining with (4.7.3) yields 


D' = ———m, 8 m.. (4.7.5) 


Here, (2 denotes the constant partial molar volume associated with formation of Li-metal 
due to electrodeposition and may be calculated from experiments or ab-initio simulations. 
Given (4.7.5) and recalling that D" = € N' from (4.3.13), the direction of electrodeposition 


induced deformations is specialized as 


N" m, 89 m,. (4.7.6) 


4.7.2 Free Energy 


Next, we consider a separable free energy per unit volume of the form 


A 


Yr(C*, C, Es V€, dr, d, Vd) = Pr (C°, Ç; d) + Wa (C) a YEE) + Wes (|VE|) a5 


+ we(C®, E, de) + 1e(|Vdl) 


(4.7.7) 


All individual functions in (4.7.7) will be specialized as isotropic functions of their argu- 


ments. Here: 


(i) YE (C*, €,d) is the contribution to changes in the free energy density due to elastic 
deformation of the host material. A consequence of isotropy is that the free energy 


wr may be expressed as (cf. Anand [113]), 


Pe (C°, £, d) = Yr (Ej, Ez, Es, E, d), (4.7.8) 


where (Ef, ES, E$} denote the principal values of the logarithmic elastic strain de- 
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fined as 
3 


3 
E° = nU’ = Y Ej 8 rf = X Ind rf 8 rf. (4.7.9) 
i=1 


j= 
Further, given (4.7.9) above, straightforward calculations (c.f. Anand and Su [196]) 


show that the mechanical contribution to Mandel stress" is then given by 


3 
m( pre He pe 
me = pas él E A (4.7.10) 


i=1 


We further decompose our free energy formulation due to elastic deformations into 
a “positive” part, Wm" due to tension and a “negative” part, Wm due to compression 
(cf. Miehe et al. [197], Klinsmann et al. [51], and Navidtehrani et al. [198]), which 
yields 


Ue (E°, E, d) = g(d) yz" (ES, E) + Yr (E°, E) (4.7.11) 


where g(d) is a degradation function of the damage variable, d. 


The positive” and "negative” mechanical parts of the free energy are individually 
specialized to the classical strain energy function of isotropic linear elasticity to 
moderately large elastic deformations using the logarithmic strain measure, E° (c.f. 


Anand [199]). They are given by 


ye) = J 


ee (ent ae en?) + GG 90) G +E + es) | 


vm (BY, £) E Je 


ale + +92) +5 (x0-300) (E+E +.) 
(4.7.12) 


where, G(£) and K(£) denote the reaction-dependent shear and bulk moduli respec- 


tively. Consistent with (4.7.9), {E;|i = 1, 2,3} denote the principal elastic logarith- 


*We refer to this as the “mechanical” contribution to the Mandel stress as it is based only on the derivative 
of the mechanical free energy, Yg’, and not the entire free energy, Wz. As such, it omits certain possible stress 
contributions, such as the Maxwell stress due to electrostatics. A detailed discussion on the significance of 
these contributions is presented in Remark 4 below. 
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mic strains. For notational convenience in (4.7.12), we additionally introduce 


ES ifE2>0 ES ifES<0 
(Ed = and (E;)_ = (4.7.13) 


O otherwise O otherwise 


We note here that only the “positive” part of the free energy 1s degraded in (4.7.11), 
while the “negative” part of the free energy remains undegraded. This in turn prevents 
damage from occurring under purely compressive stresses. The split in free energy 
due to elastic deformations is necessary not just to maintain resistance of the material 
under compression, but to additionally enable damaged regions, which have become 
electroplated (i.e. regions with d = 1 and € = 1) to develop and sustain compressive 


stresses. 


The degradation function is taken to be a monotonically decreasing function of d, 
which concurrently satisfies g(0) = 1 and g(1) = 0. Consistent with Miehe et al. 
[197], we adopt 

g(d) = (1 — d}? + €, (4.7.14) 


where e © 0 is a small positive-valued constant introduced to prevent ill-conditioning 


of the model when d = 1. Finally, the reaction-dependent shear and bulk moduli 


G(€) and K(£) obey 


G(E) = (1—h(E))C% + h(E)G%, 
(4.7.15) 
K(E) = (1—h(E))K™ + h(E)K™, 


which, consistent with (4.3.10), are interpolated using the function h(€). The specific 


form of h(£) is given by 


h(E) = E (6£? — 15€ + 10). (4.7.16) 
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This form of the interpolation function (4.7.16) is widely used in the phase-field 
literature (c.f. Chen et al. [169] and Guyer et al. [200]). Throughout, we invoke 
the “SE” and “M” superscripts to denote quantities associated with the solid-state 


electrolyte and the electrodeposited Li-metal respectively. 


(ii) ~<(c) is the change in chemical free energy due to mixing of ionic species with the 
solid host. As a simple continuum approximation to mixing, we take the chemical 


free energy to be given by a regular solution model (cf. DeHoff [201]) as 
bee) = poc + RY Cinax (e Inz+ (1-3) In(1 — 2) (4.7.17) 


Here, € = C/Cmax With € € [0, 1] represents the normalized species concentration in the 
solid host, while Cmax denotes the maximum species concentration in moles per unit 
reference volume when all the accommodating sites in the host material are filled. 
Additionally, 4o denotes the reference chemical potential of the diffusing species, R 


the universal gas constant, and Y the absolute temperature. 


(11) WEE ) includes a double-well energy function associated with the energetic barrier 
across the phases in addition to an electrochemical energetic contribution associated 
with the standard (reference) chemical state and electrostatic potential of the elec- 


trodeposited Li-metal. It is given by 
PEE) = WEC = EY + Eng + EF do. (4.7.18) 


Consistent with the works of Guyer et al. [200] and Cogswell [202], the first term 
governed by W, sets the height of the energy barrier across the phases in our contin- 
uum kinetics formulation. The second term defines the reference chemical potential, 
us of the electrodeposited species. Finally, the third term represents the energetic 


contribution associated with the electrostatic potential of the metal, o (cf. Chen et 
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al. [169] and Bucci et al. [203]). 


(iv) ve “(| VE|) is the interfacial free energy, which penalizes sharp interfaces and is sim- 
ply taken as 
1 
a (1VEl) = ¿Al VE, (4.7.19) 


with A¿ > 0, a gradient energy coefficient. This energetic contribution regularizes 
the problem from a numerical standpoint by defining a finite interface width — con- 


trolled by A¿ — for the phase-field reaction coordinate, E. 


(v) UE(CS, £, dx) is the electrostatic energy, which is taken as (cf. Narayan et al. [177] 


and Suo [204]) 
p =J (z - a) ; (4.7.20) 
Le 


Here, <(£) denotes the reaction-dependent effective electrical permittivity, and d is 
the electric displacement in the current configuration. This function is isotropic in 


the deformed body as it depends only on the magnitude of d. 


Using d = J~!Fd,, we may express (4.7.20) in terms of the referential electric 
displacement as 


1 
UC", de, E) = J (za ; ca,) (4.7.21) 
The reaction-dependent effective electrical permittivity is expressed as, 


e(€) = er(E)éo (4.7.22) 


where £o = 8.85 - 10~'* F/m denotes the electrical permittivity of vacuum and e,.(€) 
denotes the reaction-dependent relative electrical permittivity of each material phase. 


Consistent with other interpolations in this theory, we assume 


e,(é) = (1 z MEES + nee, (4.7.23) 


(vi) w4(|VdJ) is the free energy contribution associated with gradient effects of damage. 
We specialize its form as a quadratic function of the magnitude of the gradients in 
damage, |Vd| such that 

yi Va) = ¿42 |vap. (4.724) 


Here, 4* > 0 is a coefficient with units of energy per cubic volume, while £ > 0 
denotes an internal length scale, which sets the width of the damage zone across 


which the damage field varies in a diffuse fashion. 


Finally, combining (4.7.17) through (4.7.24), the total free energy per unit volume ac- 
counting for the combined effects of mixing, reaction, electrostatics, finite elastic deforma- 


tions and damage is given by 


Ye (ES, C, É, V£, dr, d, Vd) = g(d) (ES, £) T pe (ES, E) 


+ oc + RV Cmax (EINE + (1 — e) In(1 — c)) 
E as a 1 7 (4.7.25) 
+ WE? (1 — £)? + Eng + EF do + ¿AVE 
1 1 
+J! (70 . ca, + 4" [Vdl* 
Le 2 
with VP+ (E°, £) and vi (ES, £) defined in (4.7.12). 
4.7.3 Stress 
Using (4.7.25) in (4.7.10), the Mandel stress is given as 
M° = g(d)M” + M” (4.7.26) 
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with 


M“ = y lacie. (x0-360) (+B +E.) nor: 


: (4.7.27) 
= a a 
Me => [2C OEN- + (KO - FCE) (E + BS +E) [rf ort 
i=1 
and the Cauchy and Piola stress measures related by 
T=J RMR" and T, = J (R'M'R")E. (4.7.28) 


Remark 4: We note here that in deriving the stress relations (4.7.26) and (4.7.27), 
we omit the contribution to stress due to electrostatics (i.e. the Maxwell stress). That 
is, we neglect the derivative of the elecrostatic energy y£¿(C*, €, dr) with respect to 
the elastic Cauchy-Green tensor. A detailed derivation of the Maxwell stress due to 
electrostatics, leading to the classical form T* = «(e & e — (1/2)(e- e)1), with e the 
spatial electric field, may be found in the work of Narayan et al. [177]. For the case of 
inorganic solid electrolytes, as modeled here, a relative permittivity value e, = 50 has 
been experimentally reported [205]. In light of (4.7.22), the permittivity of the solid 
host evaluates then to a negligibly small number. As detailed in the works of Natsiavas 
et al. [82], Shishvan et al. [86] and Narayan et al. [177], the Maxwell stresses due to 
electrostatics are negligibly small compared to the stresses generated by the deposition 
of metal filaments inside the solid host, and can be neglected without compromising 
accuracy. As such, we omit the contribution to mechanics due to electrostatics so as not 


to obfuscate the more relevant physics discussed next. 
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4.7.4 Electrochemical Potential. Flux 


Using (4.7.25) in (4.6.6)1, the electrochemical potential of the ionic species transporting 


across the solid host is given by 


ue = po + Róln (+) + Fé. (4.7.29) 


We specialize the flux of ionic species across the solid host to be isotropic and write 
(4.6.6) as 
je = —m(c, E) VLE, (4.7.30) 


with the effective mobility, m given by 


m(e,£) = mo(E)e(1 2), and mol) = (1 — h(E) )m§ HEM 4730 


Here, m3" and mY are related to the diffusivities of each material phase through the stan- 


dard relations, 
SE 
SE _ Di 


mo” = py and my = = (4.7.32) 


with DSE and DM denoting the diffusivity of ionic species in the solid host and the elec- 


trodeposited Li-metal. 


4.7.5 Reaction Driving Force. Electrodeposition Kinetics 


Using (4.7.25) in (4.6.8)2, the electrochemical potential of the deposited Li-metal is given 
by 


= nh + (wea a?) b Foo — Jh(E)M*:N' — Div (Ae V£). (4.7.33) 


We note that in deriving (4.7.33), we neglect the derivatives of J° and the reaction-dependent 


shear, G(£) and bulk modulus K (£), present in the mechanical free energy, 4”, with re- 
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spect to the reaction coordinate, £. These terms are quadratic in the elastic strains, and 
expected to be much smaller in magnitude than the other terms in (4.7.33) (c.f. Narayan 
and Anand [178], Di Leo et al. [42]). 

Further, using (4.7.29) and (4.7.33) in (4.6.8),, the thermodynamic driving force, F for 


electrodeposition takes the following form 


F=p — pe 
E d 7 7 
= (151) -r0 (E) + Flou— 8) + E (wea). 
Soy —¢ =— — dt 
energetic _—_— electric =a (4.7.34) 


entropic energetic barrier 


— Jh(EM*:N— Div (Ag VE). 
OTN 


mechanical numerical 
regularization 


Here, the first term denotes the difference in reference chemical potentials between species 
transporting across the solid host and species in the electrodeposited metallic compound. 
The second term captures the role of configurational entropy. Critically, enrichment in con- 
centration of ionic species at the reaction sites favors electrodeposition of Li-metal, while 
depletion in concentration of Li-ions accordingly retards electrodeposition. The third term 
serves as an electric driving force for electrodeposition, stemming from the difference in 
electric potential in the electrode and electric potential in the solid conductor at the reaction 
site. The fourth term introduces a local driving force to the reaction associated with the en- 
ergy barrier between the phases, which drives the reaction towards the two energy minima. 
The fifth term captures the effect of mechanical stress on the reaction driving force. In light 
of the discussion in Sect. 4.7.1, only the stress component normal to the electrode-solid 
conductor interface affects the driving force for electrodeposition (c.f. Ganser et al. [195] 
and Deshpande and McMeeking [192]). This coupling is different in form to conventional 
diffusion-deformation theories employing an isotropic chemical distortion, in which case 
the reaction driving force couples to mechanics through a pressure term (c.f. [42, 108]). 


Critically, compressive stresses retard the electrodeposition of metallic compound at the re- 
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action sites, while tensile stresses accordingly facilitate plating. The final term in (4.7.34) 
arises from the gradient phase-field nature of the theory and regularizes the interface width, 
setting a minimum width controlled by Ae. 

Having discussed the reaction driving force, F, we now specify the reaction kinetics 
equation. Consistent with electro-kinetics theory, we invoke a Butler-Volmer non-linear 


reaction kinetics formulation and evolve the extent of electrodeposition, E as follows, 


045 (1 -a)F l = 
E R (exp (= ) — exp (S) if WWE. 
¿= ? Rd Rd (4.7.35) 


0, if €=1. 


Here, a denotes a symmetry factor, representative of the fraction of surface overpotential 
promoting anodic or cathodic reaction at the electrode interface, while Ry > 0 denotes 
a positive-valued reaction constant. For elementary single-electron transfer reactions, a 
value a = 0.5 is typically employed. Note that (4.7.35) satisfies the dissipation inequality 
(4.5.12), such that -F € > 0, where electrodeposition of Li metal, id > 0 proceeds under 
a negative reaction driving force, F < 0, consistent also with electro-kinetics theory. 

We emphasize that the focus of this work is on modeling growth of metal filaments 
across a solid conductor with continuous deposition (i.e. E > 0). As a result, we do not 
model the reverse “stripping” process (i.e. E < 0) and assume perfect contact between the 
electrode and the solid conductor is maintained at all times. Models considering the role 
of electrochemical cycling on the evolution of imperfections at the Li-metal/SSE interface 
have recently been proposed (cf. Zhao et al. [87]) that complement the proposed framework 


which focuses on Li-metal filament propagation across the SSE. 


Restrictions on electrodeposition kinetics 


We now introduce two phenomenological restrictions on the electrodeposition kinetics 


formulation presented above. First, the phase-field reaction formulation presented here 
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does not explicitly account for the presence of electrons in the electrodeposition reaction 
(4.7.1), which are necessary for Li-metal plating to occur. Since the solid-state electrolyte 
1s electronically-insulating, we assume that electrons are readily available only at the Li 
metal/SSE interface. We account for this phenomenon by restricting the electrodeposition 
reaction (4.7.35) to occur only in regions where there is existing Li-metal deposits (and 
henceforth readily available electrons), modeled as regions with € > 0. Phenomenolog- 
ically, we model this behavior by scaling the reaction rate constant, Ro with a reaction- 


dependent logistic function, f,(£) of the form 


fi = nl) — 9110), with gE = (1+exp(-arE-8))) 7. 4.736) 


The logistic function (4.7.36) is guaranteed to satisfy fı(0) = 0, such that no Li-metal 
plating will occur in the solid-state electrolyte bulk where no electrons are present. The 
parameters a, and 3; in (4.7.36) respectively control the steepness of the logistic function 
and the midpoint location. 

Second, as discussed in Sect. 4.2, we model plating to occur only within damaged 
regions of the solid-state electrolyte, which provide the necessary vacant sites to accom- 
modate Li-metal deposition. The theoretical formulation thus far does not guarantee that 
€ > 0 only in regions where d > 0. To phenomenologically enforce this constraint, we 
once again scale the reaction constant, Ry by a damage-dependent logistic function of the 


form 


fld) = o(d) — ga(0), with ga(d) = (1+ exp( — a(d 80) (4.737) 


Scaling of the reaction kinetics by the logistic function f(d) shown in (4.7.37) ensures that 
electrodeposition of Li-metal inside the solid conductor is confined only within regions of 
the host electrolyte which have incurred some degree of damage. 


Given these constraints, the reaction kinetics formulation (4.7.35) is replaced with the 
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modified form 


with the logistic functions defined in (4.7.36) and (4.7.37). 

The specific logistic functions chosen are shown in Fig. 4.2. For both functions, we 
choose a; = a 2 = 90 to ensure a steep transition. For the function f,(€), which phe- 
nomenologically models the requirement of electrons being present, we choose a small 
value of the offset point 6, = 0.05. 

(a) 1.0 


0.8 p 


o O02 04 ¿08 08 1.0 


Figure 4.2: Logistic functions modulating the reaction constant Ro with (a) the function 


fi(€) and (b) the function f2(d). Reproduced with permission from [19]. 


This suffices to enforce that Li-metal plating occurs in the presence of existing Li-metal 
deposits (i.e. in the presence of electrons), without altering the reaction kinetics signifi- 
cantly. For the function f2(d), which restricts electrodeposition within damaged regions of 
the host electrolyte, we specialize G2 = 0.2. This essentially requires that ~ 20% damage is 
incurred on the solid electrolyte before Li-metal platting can occur. In Sect. 4.9 below, we 
provide a detailed discussion on the role of the logistic the function, f2(d) which couples 


reaction kinetics and damage. 
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4.7.6 Electrostatics 


Using (4.7.25) in (4.6.9), the referential electric field is given by 
Ly 
a a Cdk. (4.7.39) 


Invoking the relations e = F~‘e, and d = J~'Fdkg, one may also relate the spatial electric 


flux density to the spatial electric field as 


4.7.7 Damage 


The evolution equation for damage is adopted from the works of Anand et al. [48] and 
Miehe et al. [197]. For brevity, we omit detailed derivations and refer the reader to the 
aforementioned publications. We note that using (4.7.25) in (4.6.10) and defining the dis- 
sipative microstress Wis AS 


Wais = Y* +Td, (4.7.41) 


the final evolution equation for damage takes the form 


Td = 2(1 — d)H — y*d + y*@Ad. (4.7.42) 


Here, I’ > 0 denotes a small viscous regularization parameter, introduced to impart stability 
to the numerical solution scheme. Consistent with the derivations in Anand et al. [48], in 
the rate-independent limit (I” = 0), the energy dissipated per unit volume as damage, d 
increases from 0 to 1 evaluates to /*. Naturally then, in our formulation, 4* represents 
an energy per unit volume dissipated during the fracture process. As later discussed in 
Sect. 4.9, ~* can be related to the experimentally reported fracture energy of the material 


of choice, Ge through a characteristic length scale parameter, /, which controls the diffuse 
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fracture zone. 
Finally, H represents a monotonically increasing history field function, which ensures 


that the damage irreversibility constraint (4.2.3) holds. It is given by 


H = max [R (E°, €) — 4*/2)). (4.7.43) 


s€(0,t] 


Consistent with (4.7.11), only the mechanical energetic contribution associated with tensile 


strains, 4; *, contributes to the evolution of damage. 


4.8 Governing partial differential equations for the specialized constitutive equa- 


tions. Boundary conditions 
The final set of the governing partial differential equations consists of: 


1. Force balance, eq. (4.6.5), viz. 


Div T, + be = 0, (4.8.1) 


with the Piola stress T} given by (4.7.28)2, and by the non-inertial body force. 


2. The local mass balance for the ionic species (4.6.7), which together with the flux 


relationship (4.7.30) yields 


¢ = Div (mV) — €, (4.8.2) 


with the mobility, m, given in (4.7.31), and the electrochemical potential of ionic 


species, u°, in (4.7.29). The reaction rate, € , is governed by the PDE (4.8.3) below. 


3. The reaction kinetics, governed by (4.7.38), which in conjunction with the reaction 
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driving force (4.7.34), yield the governing PDE for the reaction coordinate 


È = Of) Ro (exp (=r) Sas (553) a 


(ug — uo) — RY (E) + F(¢y — 0) 4 q (Ea) 48.3) 


— J'h(E)M*:N' — Div (AeVE) 


Here, eq. (4.8.3) constitutes a PDE due to the presence of the Div(V£) term in 
the reaction driving force F. Following (4.7.38), the reaction kinetics also obey the 


restriction £ = 0 when € = 1, which marks the end of the reaction. 


4. Gauss’s law, viz. (4.6.9), which along with electric displacement given by (4.7.39), 


and using eg = —V@ governs the electric potential, y through 


Div (dx) = F(c—cp) with dk = —eJC'V¢. (4.8.4) 


From a computational standpoint, numerically solving for ¢ using (4.8.4) is cumber- 
some and requires resolution of thin boundary layers in the range of a few nanome- 
ters, where deviations from electroneutrality occur [188, 189]. Given the represen- 
tative dimensions of the domains of interest in this work, this makes the problem 


computationally non-tractable and expensive. 


In practice, granted the regions where concentration gradients occur are small with 
respect to the domain of interest, an assumption of electroneutrality in the system is 
often invoked in the literature to numerically solve for ¢, such that the net charge 


dk = F(c— co) = 0 (cf. [169, 172, 174, 178, 206]). 


We adopt this approach in our numerical implementation. Making use of (4.2.8), we 
consider the current density, i, to be conserved and described by a Poisson equation 


including a sink term to represent the charge loss due to electrochemical reactions 
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(c.f. [169, 172, 174]), yielding a governing equation for ¿ of the form, 
Divi, = —FE with i, = —K(E)VO (4.8.5) 


which replaces the use of (4.8.4). Here, consistent with electro-kinetics theory, 


we introduce the effective conductivity, (€), which interpolates between the con- 
ductivity of the electrode, x*E, and the conductivity of the Li metal, «™, through 


R(E) = (1— h(€)) KF + (E). 


5. Damage evolution, governed by eq. (4.7.42), viz. 
Td = 2(1 — d)H —y*d+y"2Ad (4.8.6) 


with H a monotonically increasing history function which ensures damage irre- 


versibility (i.e. d > 0) defined in (4.7.43). 


Finally, to complete the model, a set of boundary and initial conditions must be supple- 
mented. Let Sı and Sj be complementary subsurfaces on the boundary OB of the body B, 
i.e. OB = Sy US, and Sy N St, = D. Similarly, let the pairs {S3, S4}, {S5, Se}, {S7, Ss}, 
and (Sy, S10} be complementary subsurfaces as described above. For the five degrees of 
freedom (Lu, Z, p, €, d} governed by the PDEs described above, we then respectively con- 


sider the boundary conditions 


u=ú on S;, and Tnk = tk on Sa; 

¿=é on Si, and —je-ne=Jj on Sy 

d= on Ss, and —d.-n=% on Ss; (4.8.7) 
E= £ on S7, and \VE-ng=0 on Sy; 

d 


=d on Sg, and Vd:nz=0 on So. 
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Here: i) ú is a prescribed displacement and t, a prescribed traction, ii) č is a prescribed 
normalized concentration and j a prescribed flux, 111) db is a prescribed electric potential 
and cw a prescribed surface charge density, iv) € is a prescribed reaction coordinate, and v) 


disa prescribed state of damage. Consistently, the initial conditions are also prescribed as 


u(X, 0) = up(X), CX, 0) = &(X), H(X, 0) = %(X), E(X,0) = &(X), d(X, 0) = do(X) 
(4.8.8) 
The fully coupled set of PDEs (4.8.1) through (4.8.6), along with the boundary conditions 
(4.8.7) and initial conditions (4.8.8) give then an initial/boundary-value problem for the 
unknowns of displacement u(X, t), normalized concentration c(X,t), electric potential 


(X, t), reaction coordinate €(X,t) and state of damage d(X, t). 


4.9 Numerical simulations 


In this section, we detail a set of numerical simulations aimed to highlight the important 
features of our theory, namely the interplay of electro-chemo-mechanical processes that 
govern the initiation and propagation of Li-metal filaments in SSEs. While the proposed 
framework is general in nature, we specialize it towards modeling a problem of engineer- 
ing relevance to commercialization of next-generation all-solid-state batteries as detailed 
below. 

The theoretical framework specialized in Sect. 4.7 contains five degrees of freedom 
{€,é,,d,u} governed by the partial differential equations summarized in Sect. 4.8. The 
fully-coupled PDEs are solved using the finite element method in Abaqus Standard [207] 
by developing a custom user-element (UEL) subroutine following the details described in 
Chester et al. [130]. We note that implementation of the reaction kinetics PDE (4.8.3) 
is done through the so called “micromorphic” formulation, where an auxiliary variable is 
introduced to ease numerical convergence [129, 208]. The implementation is akin to that of 


Afshar and Di Leo [15] and we refer the reader to this work for further details. Importantly, 


118 


we note that use of this numerical technique does not affect the numerical results, and 
produces equivalent outputs as would be generated with a direct implementation of (4.8.3). 

We specialize our theory for a Li-metal/LLZO solid-state architecture, which has been 
widely investigated in the literature (c.f [58, 59, 60, 65, 73]). This architecture, using an 
oxide-based electrolyte, is relevant due to its superior chemo-mechanical properties (i.e. 
high stiffness and fracture toughness, high ionic conductivity, excellent stability against Li- 
metal etc.). Growth of Li-metal filaments across a LLZO solid-state electrolyte has been 
experimentally observed to occur as both a single transgranular branch, as well as through 
a network of filaments proceeding across the grain boundaries of the SSE [58, 60, 73, 74, 
169]. 

In Sect. 4.9.1, we focus on modeling the growth of a single Li-metal filament traversing 
from the anode towards the cathode end. For the set of simulations in Sect. 4.9.1, we do not 
consider the LLZO microstructure composition (i.e. the presence of grain boundaries) and 
treat the solid host as a homogenized medium. Starting with a small initial imperfection 
at the electrode interface, we demonstrate the manner in which the coupled electro-chemo- 
mechanical phenomena govern the growth of Li-metal filaments with continuous fracture 
of the SSE. Additionally, we elucidate the role that mechanical confinement (i.e. mechani- 
cal boundary conditions) plays on the rate of crack propagation versus the rate of Li-metal 
filament growth. We demonstrate the ability of our theory to qualitatively reproduce the 
experimentally observed phenomenon of crack fronts propagating ahead of Li-metal fila- 
ments, with cracks only partially filled by the electrodeposited Li-metal [20, 69]. 

In Sect. 4.9.2, we consider the role of the SSE microstructure composition by model- 
ing the presence of grain boundaries through seeding of elements with varying mechan- 
ical properties. At present, there is a lack of numerical frameworks, which can incorpo- 
rate the manner in which microstructural information affects the coupled electro-chemo- 
mechanical growth of Li-metal filaments through SSEs. We use this set of simulations to 


elucidate the manner in which damage and mechanical stresses are coupled to electrodepo- 
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sition kinetics in this theoretical framework. In particular, we investigate how damage can 
impact electrodeposition through both the thermodynamically-consistent reaction driving 
force F in (4.7.34), and the reaction kinetics restrictions in (4.7.38). Finally, we demon- 
strate the utility of our theoretical framework by modeling the manner in which the mor- 
phology of Li-metal filaments changes for a LLZO microstructure depending on the relative 


fracture strength of the LLZO grains versus the grain boundaries. 


4.9.1 Modeling growth of a single Li-metal filament across a LLZO electrolyte. Role of 


mechanical confinement 


We present here an application of the fully-coupled, electro-chemo-mechanical framework 
towards modeling the growth of Li-metal filaments across a Li7La3Zr2O12 (LLZO) elec- 
trolyte. Growth of Li-metal filaments across pre-existing defects at the metal electrode 
interface for single-crystal LLZO microstructures has been extensively reported in the lit- 
erature (c.f. [60, 73]). We demonstrate the ability of our framework to numerically repro- 
duce the phenomena of Li-metal filament growth as a single branch protruding from the 
anode towards the cathode, while fracturing the SSE in the process. Through these sim- 
ulations, we elucidate the interplay of electro-chemo-mechanical processes, which govern 
electrodeposition of Li-metal across SSEs. In particular, we elucidate the manner in which 
mechanical confinement of the simulation domain plays a critical role on the resulting frac- 
ture and electrodeposition morphology. 

We consider a two-dimensional (2D), 200 um x 200 um plane-strain simulation do- 
main, meshed with 350 x 350 quadrilateral finite elements. This domain is representative 
of a half-cell system consisting of a Li-metal anode and a LLZO solid electrolyte, see 
Fig. 4.3(a). The dimensions of the simulation domain are assigned consistent with works 
in the literature invoking a similar phase-field treatment [77, 172, 174]. To simulate the 
Li-metal anode, a 5 um thick layer (c.f. [84, 174]) is initialized in an already fully-reacted 


state at the bottom of the domain, Fig. 4.3(a). We do so by prescribing a reaction co- 
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ordinate, € = 1 across this area. To simulate the presence of a Li-filled crack-like defect 
on the electrolyte’s surface [60, 164, 165], an imperfection 8 wm in width and 2 um in 
height is seeded on the anode/SSE interface by prescribing a state of damage, d = 1 and 
reaction coordinate, £= 1, see Fig. 4.3(b). The prescribed defect dimensions are compara- 
ble to experimental works investigating the role of surface defects on metal growth across 
inorganic SSEs [60, 209]. This setup, Figs. 4.3(a,b), mimics infiltration of Li-metal in a 
pre-existing imperfection. As shown later, the Li-filled defect serves to break the symmetry 
of the model and initiate subsequent growth of Li-metal filaments across the SSE. 

A vanishing normalized cation concentration, € ~ 0 is prescribed across the anode 
layer, Fig. 4.3(c), to numerically simulate a state of complete ionic species depletion in 
the pure Li-metal phase. According to the calculations of Li and Monroe [185], there 
are 15 available sites for lithium per LazZr20,2 formula unit, 7 of which are occupied 
at equilibrium, yielding an equilibrium occupancy 0 = 1/Lmax œ% 0.5 for LizLazZr20;2. 
Consistently, we prescribe a normalized bulk cationic concentration, ¢ = 0.5 as an initial 
condition across the solid electrolyte (regions of £ = 0), see Fig. 4.3(c). The same normal- 
ized concentration, € = 0.5, is additionally prescribed as a Dirichlet boundary condition on 
the top edge of the simulation domain. As shown in Fig. 4.3(d), we set the anode potential 


at the bottom of the domain, ¢ = OV and prescribe a potential, 4 = 0.20 V on the top of 


(a) (b) | (c) (a) 


= n — 
0.0 0.25 0.5 : 0.0 0.1 0.2 


Figure 4.3: Initial configuration of the Li-electrode/solid electrolyte simulation containing 
a pre-existing defect. Contours of (a) reaction coordinate £, (b) damage d, (c) normalized 
concentration ¢, and (d) electric potential ¢ (V). Reproduced with permission from [19]. 
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the domain, consistent with works in the literature [77, 175, 176]. The prescribed electric 
potential on the top edge of the domain, i.e. the electrolyte surface, serves to drive the 
system out of equilibrium, enabling for electrodeposition of Li-metal at the electrode-SSE 
interface. Altogether, Fig. 4.3(a-d) constitute the initial, stress-free configuration of the 
Li-metal electrode-solid electrolyte assembly. 

Mechanically, roller boundary conditions are applied at the top edge of the simulation 
domain, while the bottom edge is pinned. We consider two different boundary conditions 
for the left and right edges of the simulation domain to illustrate the role of mechani- 
cal confinement on electrodeposition kinetics and resulting fracture/electrodeposition mor- 
phology. It has been experimentally reported via in-situ X-ray tomography the propensity 
for Li-metal filaments to sometimes only partially fill cracked regions of SSEs (c.f. [20, 
69]). The cited experimental works have demonstrated the potential for electrodeposition- 
induced cracks to traverse the entire SSE ahead of the Li-deposits, which trail behind. To 
numerically reproduce this phenomenon, we apply first a roller boundary condition to the 
sides. Subsequently, we relax this lateral mechanical constraint and consider the sides to 


be free. 


Remark 5: Note that the region being modeled as the Li-metal anode shown in 
Figs. 4.3(a,b) is characterized by a state of £ = 1 and d = 0, unlike the lithium-filled 
defect, which is characterized by € = 1 and d = 1. As noted in Sect. 4.2, in develop- 
ing the theoretical framework, we restrict Li-electrodeposition inside the SSE to occur 
only within damaged regions (i.e. £ > 0 only if d > 0). Since the Li-metal anode is 
initialized before the start of the physical simulation (i.e. in a stress-free configuration), 
we do not require that d > 0 in this region. The Li-metal anode region is characterized 
by € = 1 and d = 0 so as to model a solid Li-metal layer capable of sustaining both 
tensile and compressive loads, as would be physically present in a Li-metal/LLZO half 
cell. In addition, we note here that we allow for electrodeposition of Li-metal at the flat 


horizontal anode/SSE interface, which as discussed above are characterized by a state 
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of € = 1,d = 0. As such, in those regions, we omit the logistic function f2(d) in the 


reaction kinetics (4.7.38), such that electroplating on the flat horizontal interfaces may 


occur. 


The electro-chemo-mechanical properties for a Li-Metal/LLZO solid-state system are 
tabulated in Table 4.2. We emphasize that the physics based thermodynamically-consistent 


Table 4.2: Material properties for modeling of metal filament growth in a Li-Li,LazZr30 > 
system. 


Parameter Value Source 
Electro-Chemical  KkiLzo 4.43-107? S/m Li ef al. [210] 
Kli 1.0-107 S/m Chen et al. [169] 
Dt1z0 1.0-107 2 m?/s Tian et al. [77] 
Dii 1.0-10715 m?/s Tian et al. [77] 
Gis 4.22-10* mol/m? Tian et al. [77] 
Reaction Kinetics Ro 0.1 s7} Yuan et al. [206] 
Ori 1.3-10-° m*/mol  Tantratian et al. [174] 
One 0.3 This Work 
W 1.18-10% J/m? Tian et al. [77] 
Ae 8-10* Numerical 
pJ-jm°/pmol? 
Mechanical Eizo 150 GPa Narayan and Anand [84] 
VLLZO 0.26 Narayan and Anand [84] 
Eji 4.91 GPa Narayan and Anand [211] 
Vij 0.36 Narayan and Anand [211] 
Damage T 0.04 MPa-s This Work 
yr -£ 6.67 J/m? Narayan and Anand [84] 
£ 5.0 um Regularization Parameter 


formulation allows virtually all material parameters to be found from the literature, either 


experimentally or from ab-initio simulations. That is, the theoretical framework may be 


123 


calibrated in a straightforward fashion and as such should prove useful in application to a 
number of engineering problems of relevance. 

As detailed in Sect. 4.7.2, the parameter, 7)" represents the dissipated energy per unit 
volume during the fracture process, while / denotes a length scale to account for the damage 
gradient effects. As discussed by Narayan and Anand [212] and De Borst and Verhoosel 
[213], the gradient damage formulation will be mesh-independent, provided the element 
size, he, is sufficiently small compared to the length scale, ¢. Typically, an element size 
he < 0.2 - £ is considered sufficient. Additionally, the product 4* - Z is related to the 
macroscopic critical energy release rate, Ge through Ge = Yy* -  [48, 55, 56]. Provided 
the length scale parameter, £ is chosen in a suitable physically realistic range, ¢ acts as an 
adjustable regularization parameter in our gradient damage theory. We prescribe a length 
scale, £ ~ 5 um to obtain a width for the fracture zone comparable in dimensions to pub- 
lished works in the literature [174, 206]. The characteristic element size h. = 570 nm is 
then chosen and satisfies the condition he < 0.2 - £ for mesh-independent results. This 
element size dictates the choice of a 350 x 350 finite-element mesh for the given 200 um 
x 200 um simulation domain. This ensures the interface is sufficiently narrow, while the 
number of elements required to resolve the interface does not become computationally in- 
tractable. With £ fixed, the value of /* is then chosen such that 4* - 1 = G. = 6.67 J/m? 
for LLZO. 

Additionally, we note here that choice of the viscous regularization parameter, I’ is done 
so that it has a negligible effect on the simulation results across the various cases considered 
in this work, while imparting numerical stability to our solution scheme. The value of the 
viscous regularization parameter for the framework at hand is also consistent with similar 


works in the literature (c.f. Narayan and Anand [212] and Miehe et al. [214]). 


Remark 6: To mitigate the generation of large stresses due to the elastic (rather than 
elastic-viscoplastic) model employed here for Li-metal, we reduce the overall reaction- 


induced expansion associated with deposition of Li-metal within cracks. We assume a 
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30% expansion of Li-metal during electrodeposition in the m, direction. That is, we 
choose Qué™* = 0.3 and with Qui = 1.3 - 10—*m*/mol known experimentally (cf. 
Tantratian et al. [174]), we use €™* = 2.31 - 10% mol/m?® in our simulations. The- 
oretically, one should choose (2,¡£"** = 1 (cf. Chen et al. [169]). We purposely 
use a smaller total expansion during electrodeposition to limit stress generation in Li- 
filaments. The value of 0,¡£"*% = 0.3 was calibrated to reproduce stresses within 


Li-filaments comparable to theoretical predictions based on linear elastic fracture me- 


chanics models (cf. Klinsmann et al. [164]). 


The onset and evolution of the Li-metal filament with continuous fracture of the LLZO 
electrolyte is shown in Fig. 4.4. For completeness, we illustrate the evolution of all five 
fields in time, namely from top to bottom (a) the extent of electrodeposition €, (b) the 
normalized concentration ¢, (c) the electric potential @, (d) damage d, and (e) the horizontal 
stress component, T+. Note that the first column of Fig. 4.4 shows the initial configuration 
of the Li-metal electrode-SSE assembly, where the domain is stress-free. 

Starting with a small imperfection at the anode-electrolyte interface, we can observe the 
growth in time of the amplitude of the Li-metal filament relative to the flat electrode sides, 
Fig. 4.4(a). Fig. 4.4(b-c) illustrate the evolution in normalized concentration and electric 
potential. Critically, concentration of Li-ions remains unchanged in the bulk of the SSE, 
where the unaffected electrolyte (i.e. € = ĉo where € = 0) preserves a state of electroneu- 
trality. In contrast, significant concentration gradients arise at the Li-Metal/SSE interface, 
where Li-ions are continuously consumed to plate fresh layers of Li-metal. Owing to the 
much higher conductivity of the Li-metal electrode compared to the SSE, the electric poten- 
tial across the freshly plated regions strictly follows the anode potential. These profiles are 
consistent with works in the literature invoking a similar phase-field framework to model 
the growth of Li-metal filaments across SSEs [58, 174]. 

Growth of a dominant Li-metal filament — in contrast to uniform plating of Li-metal 


— is attributed to localization of reaction kinetics at the tip of the protrusion, a phenom- 
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ena governed by both electro-chemical and mechanical forces. From an electro-chemical 
standpoint, growth of Li-metal filaments is directly related to the concentration of Li-ions 


and electric potential at the reaction site, as described by the reaction driving force, F in 
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Figure 4.4: Evolution in time for (a) the extent of electrodepostion E , (b) the normalized 
concentration c, (c) the electric potential ¢(V), (d) damage d, and (e) the horizontal stress 
component Tzs. The simulation is mechanically fully-constrained, with rollers present on 
the left and right edges. Reproduced with permission from [19]. 
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(4.7.34). Fig. 4.4(b-c) illustrate the localization of larger concentration and electric po- 
tential gradients at the tip of the Li-metal filament, as it propagates from the metal anode 
towards the far electrolyte end. This results in an increase in driving force for electrodepo- 
sition, further localizing reaction kinetics and speeding up metal filament growth at these 
sites. 

From a mechanical standpoint, Fig. 4.4(e) illustrates the evolution of the horizontal 
stress stress, T,,, across the simulation domain. Li-metal deposition inside the crack in- 
duces a large build-up in normal stresses acting on the crack surfaces. Owing to the stiffer 
nature of the solid electrolyte, Li-metal confined on the crack interior undergoes large com- 
pressive stresses, which in turn retard electrodeposition at the crack edges through the me- 
chanically coupled reaction driving force F in (4.7.34). The reverse is true at the tip of the 
crack, where large tensile stresses develop as the crack sides expand to accommodate the 
deposition of new lithium. Tensile stresses at the crack tip facilitate electrodeposition of 
Li-metal, further driving filament formation, and importantly inducing successive damage 
of the solid-state electrolyte. As shown in Fig. 4.4(e), this state of stress is sustained as 
large tensile stresses drive crack propagation and in turn allow for further Li-metal plating 
inside the cracks, which sustains filament growth. 

Finally, Fig. 4.4(d) illustrates the propagation of the electrodeposition-induced crack 
with simulation time. We emphasize that we do not predefine the crack path in our sim- 
ulations. Instead, cracks are free to evolve in any arbitrary orientation, driven by a ther- 
modynamically consistent driving force (4.7.42). Here, as there are no asymmetries in the 
simulation domain, the crack and subsequent Li-filament grow in a straight fashion from the 
anode through the electrolyte. The numerically simulated crack morphology is consistent 
with experimental observations of crack propagation in single-crystal SSE microstructures 
[20, 69, 73]. 

The importance of mechanical stress on the crack-electrodeposition interplay can be 


demonstrated by varying the degree of mechanical confinement to the solid-state archi- 
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tecture. Fig. 4.5 shows an overlay of the isocontours of the damage field, d (black lines) 
and the reaction coordinate, € (red lines) for both a simulation with (a - top row) left and 
right sides of the domain constrained, i.e. ”laterally constrained’, and (b - bottom row) 
laterally free sides. The simulation results are shown at three synchronized times denoted 
{ti, t2, t3}. We note here that the laterally free simulation terminates at t2 when the crack 
reaches the top edge of the simulation domain, while the laterally constrained simulation 


continues. 
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Figure 4.5: Evolution in time of the isocontours of damage d (black lines) and reaction 
coordinate £ (red lines) for (a - top row) a laterally constrained simulation, and (b - bot- 
tom row) a laterally free simulation. Comparing (a) and (b) demonstrates the critical role 
of mechanical confinement in dictating the degree to which the electrodeposition-induced 
crack becomes filled with Li-metal. Reproduced with permission from [19]. 
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We can clearly observe a significant difference in the extent of the Li-metal plated 
crack when comparing the laterally constrained, Fig. 4.5(a), and laterally free, Fig. 4.5(b) 
simulations. In the laterally constrained case, the damage and extent of reaction fields 
evolve in a one-to-one fashion, with the crack remaining filled with Li-metal throughout. 
In contrast, the laterally free simulation predicts that the electrodeposition-induced cracks 
propagate significantly ahead of the growing Li-metal filament. In Fig. 4.5(b), we can 


observe the electrodeposition induced crack traversing the entire solid electrolyte domain, 
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while the Li-metal filament trails behind. These results are in good qualitative agreement 
with the experimental observations of Ning et al. [20] and Hao et al. [69], where in-situ X- 
ray tomography revealed the propensity for Li-metal filaments to sometimes only partially 
fill cracked regions of the SSE, with cracks traveling at a much faster rate to reach the 
opposite cathode end plenty in advance of Li-deposits. 

Next, in Sect. 4.9.2, we demonstrate the manner in which microstructural features 
may be incorporated into our framework. Through this modeling, we elucidate the role 
of mechanics and damage on electrodeposition kinetics and resulting morphology of Li- 
metal filaments. Additionally, we demonstrate the utility of the framework by modeling the 
manner in which the morphology of Li-metal filaments changes for a LLZO microstructure 


depending on the relative fracture strength of the LLZO grains versus the grain boundaries. 


4.9.2 Modeling growth of Li-metal filaments in the presence of microstructure heterogeneity. 


We present here a set of numerical simulations 


aimed to illustrate the manner in which grain 
boundaries may be modeled within the proposed 


framework. Through these simulations, we eluci- 


date the role of mechanics and damage on elec- boundary sites, Yg, 
Initial Li-Metal 


trodeposition kinetics and resulting morphology of 


Figure 4.6: Schematic of the simu- 
lation domain highlighting elements 
electrolyte. In Sect. 4.9.1, we considered only a chosen to represent grain boundaries 
in the LLZO solid-state electrolyte. 
Reproduced with permission from 
[19]. 


Li-metal filaments as they grow across a LLZO 


homogeneous LLZO electrolyte, naturally result- 
ing in the growth of a straight Li-metal filament 
through the SSE. 

We consider here a simple representation of a polycrystalline LLZO microstructure as 
schematically illustrated in Fig. 4.6. The hexagonal pattern illustrated in Fig. 4.6 shows the 


elements within the simulation domain, which represent the grain boundaries. A character- 
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istic grain size of 30 um is chosen for the design, consistent with the experimental reports 
of Sato et al. [215] and Sharafi et al. [216]. These grain boundary elements represent 
weaker structural sites characterized by a lower fracture energy, y, compared to the bulk 
SSE, which is characterized by a fracture energy, fx As such, grain boundaries are sus- 
ceptible to damage and subsequent Li-metal electrodeposition [58]. All other mechanical 
properties are maintained the same for the entire simulation domain. 

From an electro-chemical standpoint, we assume the grain boundaries to have the same 
electro-chemical properties (i.e. diffusivity, conductivity) as the grain bulk. This assump- 
tion is purposely chosen to enable us to investigate purely the mechanical effects of grain 
boundaries on electrodeposition kinetics and resulting Li-metal filaments morphology. As 
such, any effects sourcing from variation in electro-chemical properties between grain 
boundaries and the solid electrolyte bulk are neglected, although we note that these effects 
could also be incorporated in the theoretical framework presented. All material properties 
remain the same as those in Sect. 4.9.1 and listed in Table 3.1. 

Boundary conditions are identical to the fully-constrained simulation described earlier 
in Sect. 4.9.1 and consistently here, an electrical potential, 4 = 0.2 V is applied to drive 
the system out of equilibrium. We employ a similar configuration for the Li anode-SSE 
interface and introduce a small Li-filled defect. With Fig. 4.6 in mind, we can observe that 
as crack and Li-metal filament grow across the SSE, they will eventually intersect with the 
mechanically weaker grain boundaries. At this point, Li-metal filaments can either propa- 
gate: 1) intergranularly across the mechanically weaker grain boundaries, ii) transgranularly 
by fracturing the bulk LLZO and continuing to grow in a straight line, or iii) through a com- 
bination of both modes. The growth direction and resulting morphology of Li-filaments is 
dictated by the thermodynamic reaction driving force (4.7.34), which models the role of 
chemical, electrical and mechanical driving forces on Li-metal electrodeposition kinetics. 

To illustrate the manner in which mechanics couples to electrodeposition kinetics, we 


perform a series of numerical simulations for the specific case of Yj, = (1/3) - Wu 
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To distinguish the two ways mechanics couples to electrodeposition, we consider first the 
reaction kinetics (4.7.38), which may be written as 3 = fold) £,(€, F). This highlights 
the first mechanism being the logistic function, f2(d) which directly couples the reaction 
kinetics to mechanical damage. As detailed in Sect. 4.7.5, the function, f2(d) is introduced 
phenomenologically to restrict electrodeposition strictly within damaged regions across 
the SSE. Second, the electrodeposition driving force F from (4.7.34) may be written as 
F = Fs(c, E, ¢) — J“h(€)M* :N', where we now separate the stress-dependent term, which 
couples the reaction driving force to mechanical stresses. 

Fig. 4.7(a) illustrates the fully-coupled simulation, where both mechanisms are active. 
We can observe the role of mechanically weaker grain boundaries in guiding the crack and 
subsequent growth of Li-metal filaments across the SSE. As shown in Fig. 4.7(a), Li-metal 
filaments growth remains almost entirely confined within the damaged grain-boundaries, 
with a small vertical crack within the first grain appearing, but failing to continue to grow. 


(a) Fully-Coupled : (b)  Stress-Coupled Only 
F = Fry(c,€, 6) — JRM: N" | F = Fs(c,E, $) — JUh(EM*:N' 


(c) Kinetics-Coupled Only 
E= fo(d)Fi(E, F) 
F= Ê, (c, E $) 


Reaction Coordinate, € 


Damage, d 


Figure 4.7: Contours of reaction coordinate E (top row) and damage d (bottom row) for (a) 
a fully-coupled simulation, (b) a stress-coupled only simulation, and (c) a kinetics-coupled 
only simulation. Reproduced with permission from [19]. 
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Fig. 4.7(b) shows a numerical simulation in which the logistic function, f2(d) is re- 
moved, and direct coupling of electrodeposition kinetics to damage is suppressed. The 
only coupling between mechanics and electrodeposition is through the stress-coupled reac- 
tion driving force F. Importantly, we can observe that even in the absence of a restriction 
preventing plating outside of cracks, Li-metal filaments propagate primarily in cracked re- 
gions due solely to the stress-state within the cracks creating a favorable reaction driving 
force. Consistent with the results in Fig. 4.7(a), we observe the tendency for Li-metal to 
preferentially deposit across the mechanically weaker grain boundaries, even though a con- 
current transgranular branch forms in this case. This result illustrates that the presence of 
the logistic function, f2(d) in (4.7.38) does not significantly alter the resultant morphology 
of Li-metal filaments across the SSE microstructure. 

Finally, Fig. 4.7(c) illustrates the reverse coupling, where the stress-coupling of the 
reaction driving force F is suppressed, while the reaction kinetics remain coupled to dam- 
age through the function f2(d). In Fig. 4.7(c), we observe an entirely different pattern of 
crack propagation and Li-metal deposition across the SSE. In contrast to simulation results 
with active stress-coupling, Figs. 4.7(a,b), in which cracks and Li-metal filaments predom- 
inantly evolve across the weaker grain boundaries, the simulation with no stress-coupling 
predicts cracks and Li-metal filaments to grow primarily in a transgranular fashion. That 
is, in the absence of a stress-coupled reaction driving force F, the presence of weak grain 
boundaries is virtually unseen by the electrodeposition kinetics. This difference in pattern 
is expected, granted in the absence of stress-coupling, electrodeposition kinetics is dictated 
by purely electrochemical forces. As discussed in Sect. 4.9.1, localization of higher concen- 
tration and electric potential gradients at the tip of the straight Li-metal branch compared 
to the slanted grain boundaries favors reaction and speeds up the growth of a dominant 
vertical Li-metal filament across the SSE. 

These results highlight the importance of coupling the electrodeposition driving force, 


F to mechanical stress. Importantly, mechanical stresses in the presence of microstructural 
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heterogeneities (1.e. grain boundaries, pores etc.) can significantly alter the morphological 
evolution of cracks and Li-metal filaments across the SSE. Such coupling has also been 
shown experimentally [58, 59, 217] and is reproduced here through the proposed theoretical 
framework. 

Having discussed the manner in which mechanics couples to electrodeposition kinetics, 
we now turn our attention to the role of microstructure composition, specifically the role 
of grain boundaries with varying fracture energy. We consider a range of grain boundary 
fracture energies, +, relative to the bulk solid electrolyte fracture energy, Vx» to explore 
the manner in which microstructural features with varying mechanical properties alter the 
growth morphology of cracks and Li-metal filaments across the SSE. 

Fig. 4.8 shows contours of reaction coordinate (top row) and damage (bottom row) for 
simulations with (a) 4%, = (1/5) Varo (0) Ye, = (1/3) Vur and vs, = (1/2) -Yax Note 
that all these simulations are fully-coupled, and thus the result in Fig. 4.8(b) is identical to 
that in Fig. 4.7(a). 

For both cases with 4%, = (1/5) - Yiu and Y% = (1/3) - Wur Figs. 4.8(a,b), we can 
observe the role of weaker grain boundaries in guiding Li-metal to preferentially deposit 
at these sites. Owing to their lower fracture energy, grain boundaries are easier to dam- 
age compared to the bulk SSE. This creates a favorable stress-state within the fractured 
grain boundaries, making it favorable for growth of Li-metal filaments to follow the grain 
boundary network, suppressing the growth of a vertical transgranular crack and Li-metal 
branch. These results are consistent with experimental observations, which report Li-metal 
filaments to preferentially grow along the grain boundaries (i.e. integranularly) in a poly- 
crystalline LLZO electrolyte [58, 59]. 

With increase in fracture energy of the grain boundaries relative to the bulk SSE, 

o = (1/2) - Yiu in Fig. 4.8(c), we observe a change in the resulting crack and Li-metal 
filament morphology. We observe now the formation of a distinct vertical transgranular 


branch propagating through the SSE in advance of concurrent Li-metal growth proceeding 
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Reaction Coordinate, € 


Damage, d 


Figure 4.8: Contours of reaction coordinate E (top row) and damage d (bottom row) with 
varying grain boundary fracture energies, respectively (a) V%, = (1/5) - Vino (b) V% = 
(1/3) - Yu and (c) Y, = (1/2) - Yiu Reproduced with permission from [19]. 

across the grain boundaries. As the fracture energy of the grain boundaries increases closer 
to that of the electrolyte bulk, fracturing these sites becomes more difficult, and transgranu- 
lar growth becomes preferable. This result suggests there exists a threshold fracture energy 
ratio, where the electrochemical driving force overcomes the weak grain boundaries and 
the growth morphology of Li-metal filaments transitions from primarily intergranular to 
primarily transgranular. Clearly, the limiting case of Wg, = Viu 18 identical to the simula- 
tion results shown in Fig. 4.4 and produces a single vertical Li-filament. 

At present, a microstructure description of the process of Li-metal filament growth 
across SSEs is missing and necessary. The theoretical framework summarized and spe- 
cialized in Sects. 4.6 and 4.7 alongside numerical simulations presented here show the 
first framework capable of predicting the resulting morphology of Li-metal filaments for 
a given SSE architecture, including a transition from primarily intergranular to primarily 


transgranular filament growth. 
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4.10 Concluding Remarks 


We have formulated and numerically implemented a thermodynamically-consistent electro- 
chemo-mechanical gradient theory, which couples electrochemical reactions with mechan- 
ical deformation and damage in solids. The framework models transport of ionic species 
across the solid host under diffusion/migration mechanisms and concurrent electrochemi- 
cal reaction at cracks across the solid host, where ionic species are reduced to deposit a new 
compound. Critically, the theory captures the interplay between growth-induced fracture of 
the solid host and electrodeposition of a new material inside cracks by tracking the damage 
and extent of electrodeposition using separate phase-field variables. 

While the framework is general in nature, a specialization towards modeling of Li-metal 
filament growth across SSEs for application in all-solid-state batteries is demonstrated. In 
particular, we elucidate the coupling of electrical, chemical, and mechanical processes, 
which govern the growth of Li-metal filaments across SSEs with and without account for 
microstructure heterogeneity. Consistent with experiments, we demonstrate the ability of 
the framework to numerically reproduce both transgranular and intergranular growth of 
metal protrusions across the SSE microstructure. 


Major contributions of the proposed framework include: 


e The framework enables for modeling concurrent diffusion, reaction, deformation and 
damage in solids. Importantly, the theory extends the current literature in continuum 
mechanics to model solids undergoing electrochemical reactions, which are intri- 


cately coupled to damage of the underlying solid host. 


The interplay between electrodeposition of new material inside damaged (fractured) 
regions is modeled. Consistent with experiments, electrodeposition is restricted 
within damaged regions across the solid conductor, which supply the necessary va- 
cant space for accommodating a new material. This feature provides an extension to 


current electro-chemo-mechanical phase-field formulations in the literature, which 
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do not model growth-induced fracture of SSE due to Li-deposition and indistinguish- 


ably allow for Li-metal plating across the host. 


The thermodynamically-consistent treatment enables for derivation of a physically 
motivated reaction driving force in which different contributions of energy, config- 
uration entropy, electrostatics, mechanical deformations and damage can be readily 
identified. Particularly useful then is the fact that material properties driving the 


electrodeposition kinetics can be calibrated from the literature or experiments. 


The framework elucidates the role of mechanical confinement on the crack vs. elec- 
trodeposition morphology. Under specific mechanical boundary conditions, we nu- 
merically reproduce the presence of cracks, which are only partially filled with Li- 


metal. 


Finally, we demonstrate the role of microstructural heterogeneities, in particular grain 
boundaries, in dictating the morphology of cracks and Li-metal filaments across the 
SSE microstructure. Both transgranular and intergranular growth mechanisms are 
numerically reproduced consistent with experiments, elucidating the critical role of 
mechanics on electrodeposition kinetics and the resulting morphology of Li-metal 


filaments across the SSE. 


Growth of Li-metal filaments across SSEs constitutes a critical degradation mecha- 
nism, hindering commercialization of next-generation all-solid-state batteries. It is thus 
critical, from both an experimental and modeling perspective, to understand the interplay 
between non-uniform electrodeposition kinetics and the SSE microstructure composition. 
In particular, the manner in which size and distribution of defects and microstructural het- 
erogeneities ultimately govern the growth of Li-metal filaments. The proposed framework, 
in conjunction with experiments, enhances our current understanding of the electro-chemo- 
mechanical processes, which govern the initiation and growth of Li-metal filaments across 


SSEs. From a design standpoint, the proposed framework provides then a utilitarian numer- 
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ical tool, which can help elucidate new strategies to mitigate growth of Li-metal filaments 


and enable for successful commercialization of all-solid-state batteries. 
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CHAPTER 5 
CONCLUSION 


5.1 Summary 


Development of high-fidelity numerical simulations for design and optimization of next- 
generation all-solid-state batteries (aSSBs) is of critical importance towards their success- 
ful commercialization to meet increasing energy demands. While an extensive body of 
experimental literature, reporting on the numerous degradation mechanisms in solid-state 
architectures exists, understanding the intricate coupling between electrochemical and me- 
chanical phenomena from a modeling perspective remains a topic of active research. In 
this thesis, we address this limitation and develop a series of numerical frameworks, which 
model the governing electro-chemo-mechanical processes coupled with fracture across the 
different constituents in aSSBs to assess their role on longevity and reliability of aSSBs. In 
particular, the developed electro-chemo-mechanical frameworks elucidate the degradation 
mechanisms in two critical constituents of aSSBS, namely (1) composite electrodes and (11) 
solid-state electrolyte due to growth of Li-metal filaments. To conclude, we summarize 
here major contributions associated with each part of the thesis: 

Part I: Modeling chemo-mechanical multi-particle interactions in composite electrodes for 


LIBs. 


e Formulated and numerically implemented a theoretical framework based on use of 
chemo-mechanical surface elements for modeling the stress-coupled chemo-mechan- 
ical interactions in complex electrode designs under galvanostatic charging. In partic- 
ular, the theoretical framework allows one to capture both chemical and mechanical 
interactions between active particles in battery electrodes under galvanostatic charg- 


ing conditions. Importantly, the developed framework captures chemical interactions 
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between particles, whereby current is non-uniformly distributed between active par- 


ticles based on their local stress-coupled interface reaction kinetics. 


Mechanically, the developed surface elements allow one to capture the role of me- 
chanical interfacial damage on electrochemical performance by coupling interfacial 
reaction kinetics to mechanical damage. This in turn enables for quantifying how 
loss of mechanical integrity at the interface limits reaction kinetics and overall elec- 


trochemical performance of the electrode. 


Using the developed simulation capability, we investigated the chemo-mechanical 
behavior of novel LIB electrode designs composed of hollow a-Si nanotubes. Elec- 
trode capacity loss due to interfacial decohesion (1.e mechanical degradation) was 
quantified by performing a series of charge/discharge cycles at varying C-rates. In 
doing so, we demonstrated the manner in which the proposed theoretical framework 
can be used to capture mechanically induced electrochemical capacity loss over a 


number of different charging cycles, consistent with experiments. 


Part II: Modeling of electro-chemo-mechanical behavior coupled with damage of com- 


posite electrodes for all-solid-state batteries. 


e Leveraging the simulation capability developed in Part I, chemo-mechanical behavior 
and electrochemical performance of a LiCoO>-Li¡yGeP2S¡2 composite cathode for 
all-solid-sate batteries was investigated. A key focus here involved understanding 
the role of microstructure on the state of mechanics and integrity of interfaces, in 
turn impacting the overall electrochemical performance. First, we studied the role of 
variations in chemo-mechanical properties of the composite electrode of choice, with 
a focus on SSE stiffness and volumetric contraction of active material upon lithiation. 
Second, we investigated the role of two critical microstructural features: 1) the size 
and size distribution of active particles, and 11) the active material volume fraction 


(i.e. packing density). 
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e Variations in SSE stiffness and active material contraction during lithiation were 
demonstrated to significantly impact the electrochemical performance of SSB cath- 
odes, limiting attainable capacity in the presence of mechanical interfacial delamina- 
tion. From a design standpoint, the developed theoretical framework thus provides a 
utilitarian tool for future exploration of manufacturing of composite electrodes with 
different active material-SSE pairs to enhance mechanical integrity and concurrently 


electrochemical performance. 


With respect to microstructural composition, particle size distribution was demon- 
strated to have a negligible effect on the stress build-up and potential for delamina- 
tion at the active particle-SSE interface. This finding was demonstrated both through 
simple single-particle RVEs and through complex multi-particle RVEs with different 


particle size statistics. 


Importantly, active material volume fraction, commensurate with packing density, 
was shown to dictate the build-up of interfacial stress at the active particles surface, in 
turn controlling the onset and evolution of mechanical interfacial degradation. While 
increasing the active material volume fraction naturally increases the capacity of the 
LCO-LGPS cathode under consideration, these gains can be quickly undermined by 
the exacerbation in state of interfacial delamination experienced across the electrode 
microstructure. The proposed theoretical framework can thus serve as a useful tool 


towards design of complex novel SSB cathode architectures for future applications. 


Part III: A novel electro-chemo-mechanical framework coupled with fracture for mod- 


eling of Li-filament growth in solid-state batteries. 


e Formulated a thermodynamically-consistent electro-chemo-mechanical gradient the- 
ory, which couples electrochemical reactions with mechanical deformation and dam- 
age in solids. The framework models transport of ionic species across the solid host 


under diffusion/migration mechanisms and concurrent electrochemical reaction at 
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cracks across the solid host, where ionic species are reduced to deposit a new com- 
pound. Critically, the theory captures the interplay between reaction-induced fracture 
of the solid host and electrodeposition of a new material inside cracks by tracking the 


damage and extent of electrodeposition using separate phase-field variables. 


The thermodynamically-consistent treatment enables for derivation of a physically 
motivated reaction driving force in which different contributions of energy, config- 
uration entropy, electrostatics, mechanical deformations and damage can be readily 
identified. Particularly useful then is the fact that material properties driving the 


electrodeposition kinetics can be calibrated from the literature or experiments. 


The framework elucidates the role of mechanical confinement on the crack vs. elec- 
trodeposition morphology. Under specific mechanical boundary conditions, the ca- 
pacity of the framework to numerically reproduce propagation of cracks partially 


filled with Li-metal was demonstrated, consistent with experiments. 


Finally, we demonstrated the role of microstructural heterogeneities, in particular 
grain boundaries, in dictating the morphology of cracks and Li-metal filaments across 
the SSE microstructure. Both transgranular and intergranular growth mechanisms 
were numerically reproduced consistent with experiments, elucidating the critical 
role of mechanics on electrodeposition kinetics and the resulting morphology of Li- 


metal filaments across the SSE. 


5.2 Outlook 


While progress has been made in this thesis towards modeling of the coupled electro- 
chemo-mechanical processes in all-solid-state batteries, much remains to be done. Below, 
we list potential directions for future work which have not been addressed in this thesis: 


Outlook on modeling composite electrodes for all-solid-state batteries 
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¢ In its current form, the framework developed in chapter 2 does not account for trans- 
port of ionic species across the solid conductor to reach the active material-SSE in- 
terface. Instead, an underlying assumption of fast ionic transport across the solid 
electrolyte is invoked, translated numerically to a uniform constant chemical poten- 
tial across the nodes at the interface in contact with the solid-electrolyte matrix. With 
increase in active material (i.e. packing density), aside from an increase in the role of 
mechanics on integrity of active particle-SSE interface, additional complexities such 
as blockage of ionic transport pathways arise. Blockage of ionic transport pathways 
can in turn potentially limit utilization of active material available for reaction, lo- 
calizing reaction kinetics near the cathode-separator interface. Thus, a fundamental 
understanding of the coupled role of electrode microstructural stochasticity and the 
spatio-temporal electrochemical-transport dynamics on the active material utilization 


is required. 


Owing to the stiff nature of all-solid-state batteries, mechanical integrity of the inter- 
face constitutes merely one of three critical mechanisms through which fracture lim- 
its performance in composite electrodes. Equally important here is development of an 
appropriate numerical framework which enables for modeling fracture at the i) bulk 
of active particles as they undergo swelling/contraction with Li-insertion/extraction 
and 1i) bulk of the solid conductor in response to any chemically induced deforma- 
tions associated with transport of ionic species or due to deformation of active par- 
ticles embedded in it. In developing these models, consideration must accordingly 
be paid to the underlying microstructure of both the active material and solid ionic 
conductor to numerically resolve the presence of grain boundaries (in the case of 
polycrystalline active materials or SSE), inherently existing pores etc. At present, 
a unified cohesive phase-field formulation which concurrently accounts for dam- 
age mechanisms at the active material-SSE interface, active particle bulk as well 


as bulk of the solid ionic conductor is missing. Development of a unified framework 
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1s critical and would enable for elucidating which, if any of the three underlying 
damage mechanisms, critically dictates electrochemical performance of composite 


electrodes. 


Towards design optimization of composite electrodes for improved performance, nu- 
merical simulation of high-resolution experimentally resolved composite microstruc- 
tures is critical. Nevertheless, this aspect significantly increases computational re- 
quirements and time. Towards this goal, development of deep-learning—AI models, 
which can connect continuum simulations into a physics-based data-driven multi- 
scale model are of particular interest and have attracted significant attention as a tool 
for rapid material design [218]. Through the proposed deep-learning model, data 
generated in predictive continuum simulations in the form of consecutive images can 
be used to train a deep neural network, which allows then to make computation- 
ally efficient predictions. As an example, use of deep-learning paradigms based on 
processing of consecutive images as time sequences leveraging ConvLSTM convo- 
lutional networks has been demonstrated with great success in predicting 1) fracture 
patterns in crystalline solids from atomistic simulations [218] and 11) complex me- 
chanical stress states in hierarchical composites for material design [219]. Success- 
fully adopting a scalable machine-learning model could enable for fast predictions of 
microstructure evolution with electro-chemical cycling of more complex electrode 
architectures composed of hundreds of particles, while bypassing the need for com- 


plex, computationally-intensive simulations. 


Outlook on modeling of growth of Li-metal filament across solid-state electrolytes 


e Continuous deposition of Li-metal inside defects induces large compressive stresses, 
which in turn can cause Li-metal to flow plastically. While the role of plastic flow 
on propagation of Li-metal filaments across SSEs is a topic of active research, an 


elastic-viscoplastic extension of the framework is a critical ingredient in modeling 
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the morphological evolution of voids with continuous plating/stripping for nucle- 
ation of Li-metal filaments at the anode-SSE interface. Due to the viscoplastic na- 
ture of Li-metal, plastic deformation and creep become relevant at sufficiently high 
pressures and can act to close the voids, reducing contact loss at the anode-SSE inter- 
face. Reducing contact loss could then mitigate the propensity for Li-metal filaments 


growth. 


In line with the scope of the framework targeting propagation of Li-metal filaments 
across the solid-state electrolyte, at this stage we purposely specialize the framework 
to model only forward reaction, i.e. Li-metal deposition. Nevertheless, we empha- 
size here that the proposed framework, through the employed Butler-Volmer reac- 
tion kinetics and the thermodynamically-consistent electrodeposition driving force 
enables for modeling both plating and stripping processes. Numerically investigat- 
ing the effect of electrochemical cycling becomes critical during the nucleation stage 
in modeling the evolution of the Li-metal anode-SSE interface with the potential for 


void formation, and constitutes an important aspect left to future explorations. 


From a design standpoint, it is important to explore different mechanisms through 
which one can tailor the solid electrolyte microstructure composition to prevent Li- 
filament growth. In this context, a natural extension and utilitarian study leveraging 
the proposed framework includes an understanding of the role of stiff mechanical in- 
clusions as a potential mechanism for hindering Li-filament growth with fracture of 
the solid electrolyte. Key to investigations of this nature is determination of an opti- 
mal inclusion size and percent infill such that one effectively limits dendrite growth, 


without significantly limiting transport pathways. 


144 


Appendices 


APPENDIX A 
DETAILS ON THE NUMERICAL IMPLEMENTATION OF 
CHEMO-MECHANICAL SURFACE ELEMENTS 


We present here details of the numerical implementation of the constitutive model for the 
chemo-mechanical surface elements described in chapter 2. In Sect. A.1, we first introduce 
the variational formulation of the surface elements, followed by specialization for a 2D 


quadilateral element in Sect. A.2. 


A.1 Variational Formulation 


We present here details on the variational formulation of the proposed surface elements. 
Note that although our surface elements are chemo-mechanical in nature, we present here 
only the variational formulation for the chemical counterpart of the elements, responsible 
for modeling of the non-linear reaction kinetics at the interface. The variational formulation 
for the mechanical portion is standard in the literature, and the reader is referred to the 
works of Park and Paulino [220] and Camanho and Davila [127] for further details. 

We choose the chemical potential, yu as the degree of freedom governing mass balance. 
The strong form of the local species mass balance, supplemented with appropriate bound- 


ary conditions is then given by, 


ċr = —Jdivj, in B, 


Mass Balance = 4 u = ji, on S,,, (A.1.1) 


—mgradu -n =j, on S; 


where B denotes the body in the deformed configuration, while S, and S; are complemen- 


tary subsurfaces of the boundary OB of the body BG in the sense that 0B = S, (US; and 
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Su MS; =I. 

Let w be a virtual variation in chemical potential, yy such that the variation w vanishes 
on §,,. The weak form corresponding to Eqns. (A.1.1) is obtained by multiplying Eqn. 
(A.1.1), by w and integrating, which yields 


[ went) ae — f wdiv(mgradu) dv = 0 (A.1.2) 
B B 


where the “grad” and “div” operators are with respect to the current configuration, and in 
writing (A.1.2) we have used the constitutive equation j = —mgradu. Application of the 


divergence theorem to the second term in (A.1.2) leads to 


f vosivtmerasy) dv = I w(mgradu) - nda — J meran) - gradw dv. (A.1.3) 
B ƏB B 


Use of (A.1.3) in (A.1.2), gives 


[went du + f (meray - gradw dv — f w(mgradu) - nda = 0 (A.1.4) 
B B aB 


We now consider the body B to be composed of two sub-domains, namely B+ and B”, 
separated by an interface Z, as illustrated in Fig. A.1. Application of (A.1.4) individually 
to B* and B- then yields 


A wér (3) av f (mgradj1)-gradw a- f w(mgrady)-nda~ f (¡"-n*)w*da = 0 
B+ B+ aB+ T+ 
(A.1.5) 


f win) dv f (mgradu)-gradw w- | w(mgrady)-nda— f (j -n~ )w da = 0 
f Be cs 
(A.1.6) 


The unit normal, n* to the surface Z” points from B* to B~, while n7, the unit normal to 
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T- points from B~ to Bt, and for notational convenience we introduce n” = —n* = ñ. 


Equivalently, the interface flux j* represents species flux across the interface from one side 


to the other along the n~ directions. 
Making use of B = B*|]B”, 0B = OB* |JOB-, T* = T and the continuity condi- 
tion, j*-n* = —j--n” = jo: at the interface, we may combine (A.1.5) and (A.1.6) to 


yield 


[wen du + f (merady - gradw du — f w(mgradu) - nda + le [[w]] da = 0 
B B ðB T 
(A.1.7) 
where [[w]] = wt — w” denotes the jump in virtual chemical potential at the interface. 
Introducing the boundary condition (A.1.1)3, noting that w vanishes on S,,, we may 
finally write (A.1.7) as 


[ went) a -+ J merani) - gradw du +/ wj da + [ic [[w]] da=0. (A.1.8) 
B B T 


as; 


Figure A.1: Schematic illustration of the cohesive interface for numerical implementation 
of galvanostatic charging conditions. Reproduced with permission from [17]. 


The terms in (A.1.8) fall under the following categories. The first three terms in (A.1.8) 
arise from the conventional mass balance across a body B without a cohesive interface. In 


turn, the last term denotes the contribution from the presence of an interface. In simulating 
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a body with an interface such as the one shown in Fig. A.1, the first three terms in (A.1.8) 
are implemented through traditional continuum elements, which discretize the parts B+ 
and B”, while the last term is implemented through the proposed surface elements, which 
discretize the interface Z. The implementation of bulk continuum chemo-mechanical ele- 
ments is discussed in detail in Di Leo et al. [129] and Chester et al. [130]. We thus focus 
the remainder of our discussion on the last term in (A.1.8), which is implemented through 
use of surface elements. This entails thus a constitutive law between the flux of chemi- 
cal species at the interface, j.. and a jump in the chemical potential [[u]], modeled under 
Butler-Volmer interface kinetics as outlined in (A.1.14). 

From a finite-element modeling standpoint, evaluation of the consistent residual and 
tangent stiffness for the chemical counterpart of the theory is then similar to conventional 
numerical implementation of mechanical cohesive elements reported in [220] and will be 
briefly summarized here for convenience. We outline the numerical implementation in 
detail, specializing for a 2D quadilateral element in Sect. A.2. 

In a finite element framework, one can compute the vector of nodal jumps in chemical 
potential at the interface of surface elements, Ap (ie Ap, Ap - c.f Fig. A.2) from 
the vector of nodal chemical potentials, yz (c.f Fig. A.2) making use of the following 
relationship in matrix form 


Ap = Ly (A.1.9) 


Analogous to formulation of the mechanical cohesive element, here L defines a local sep- 
aration matrix. It relates the nodal jumps in chemical potential at the interface to nodal 
chemical potentials of the surface elements [220]. 

In conjunction with the shape function matrix N (c.f (A.2.5)), local jump in chemical 
potential anywhere along the cohesive interface, Ay is interpolated from the nodal chemi- 


cal potential jumps at the interface of the surface elements, Ap 


Au =NAp (A.1.10) 
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Combining (A.1.9) and (A.1.10), one obtains the relationship between the jump in chemi- 
cal potential along the cohesive interface and the nodal chemical potentials of the surface 
elements, 


Au=Bp with B=NL (A.1.11) 


Similarly here, B defines a matrix analogous to global displacement-separation relation 
matrix for the conventional mechanical cohesive element [220]. 
Element level residual, Roh, and tangent stiffness, Keon, for the surface elements can 


then be computed as follows, 


R= f B” jez da (A.1.12) 
T 


_ ORM =- [st Ojez O(Ap) da=- f T ez 
== LA Bau) ðn > An) 


where flux of ionic species at the interface, je- is computed as a function of the jump in 


coh __ 


B da (A.1.13) 


chemical potential, Ay using Butler-Volmer kinetics 


. i io. .,1Ap 
cz — =? h 
J F sin) 


(A.1.14) 


A.2 Numerical Solution Methodology 


We implement our theory for modeling of multi-particle interactions in composite elec- 
trodes and develop a chemo-mechanical surface user-element (UEL) for ABAQUS/Standard 
finite element package. This section outlines the numerical implementation procedure, spe- 
cializing on a 2D-quadilateral element with three degrees of freedom (u1, uz, u) active at 
each node. 

At the initiation of each solution step, the user has access to nodal coordinates (CO- 
ORDS) and updated global degrees of freedom (U), while all specifications regarding nodal 
connectivity (NODE) and material properties (PROPS) are prescribed via an accompany- 


ing input file. As part of the user element subroutine (UEL), the user is required to eval- 
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uate/update the right-hand side matrix (RHS) and the Jacobian matrix (AMATRX) at the 


end of each solution step for ABAQUS to access whereby 


¢ The RHS matrix contains the contributions of an element to the right-hand-side vec- 
tor of the overall system of equations. Specializing for our surface element, RHS is 


the overall elemental residual which in matrix form is given by 


R = [R}_, R1, R}, R2, R2,, RRA, Re, 


n]T 
Ri] 
with n denoting the total number of nodes pertaining to an element. 


e AMATRX contains the contribution of an element to the stiffness of the overall sys- 
tem of equations. Specializing for our formulation, AMATRX is the global stiffness 


tangent given by 


a A Meas. Wade a A nai az Maa 
win Kano Kon Kana Kuzo Kisu uzn Kuzua Kuan 
Kon Koa Ku Koa Ky Ku Kon Ku Kuy 
Kon Ka Koi Kaua Kiis Ka iii Kita Ku Ko: 
ga |Een Kimo Kiu Kina Kim Kin Ka Kimo Kim 
Kon Kua Ku Koa Ko Kys Kin Kpu Kuu 
Krt Krt Krt K”? K”? K”? Krr Krr Krr 
uiui u1u2 up uiui u1u2 up uiui u1u2 up 
a Ps ap AG a Kin Kam Kums Kuzu 
Kin Kia Ki Koa Kiu Kup Kim Ku Kyuu 


While the proposed surface element is of chemo-mechanical nature, we detail here the nu- 
merical implementation of the chemical counterpart of the theory only. For the interested 
reader, a more detailed description on the numerical implementation of a mechanical co- 
hesive element under mixed-mode softening can be found in the work of Park and Paulino 


[220]. 
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Figure A.2: (a) Surface elements integrated at the boundary of electrode particles to capture 
variations in current density along a particle surface and the associated current distribution 
across different particles. (b) Mechanism of current distribution at the electrode surface as 
a function of jump in chemical potential at the particle-electrolyte interface. Reproduced 
with permission from [17]. 


Given the nodal chemical potentials of the surface elements (41, H2, 43, 44) at the ini- 
tiation of a solution step, nodal jumps in chemical potential (A pe, Ap?) at the interface 


of the surface elements are computed using 


Ap = p-m and Au = ps — pa (A.2.1) 


In matrix notation, nodal jumps in chemical potential at the interface of the surface ele- 


ments, Á yu can be obtained from the nodal chemical potentials, pz through 
Ap = Lu (A.2.2) 


with 


1001 
L= (A.2.3) 
0 <1 4 o 


In conjunction with the shape function matrix N, jump in chemical potential, Ay anywhere 


along the interface is interpolated from the nodal jumps in chemical potential at the inter- 
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face of the surface elements, Á p using, 
Au =NAp (A.2.4) 


with 


e 1 
N, = — aid: N= is (A.2.5) 


and € denoting a normalized coordinate along the undeformed surface element. 
Combining (A.2.2) and (A.2.4), jump in chemical potential along the interface of sur- 


face elements, Ay can be related to the vector of nodal chemical potentials, jz through 


Au=Bp where B=NL= [-N, — Na Na Ni] (A.2.6) 


One can then numerically compute the element level residual, Rem and tangent , Keon 
from Eqns. (A.1.12) and (A.1.13), employing a Gauss-quadrature scheme to evaluate the 
integral expressions. In our implementation, we use a two-point Gauss-quadrature rule 


with location of Gauss integration points at € = +1/3. 


A.3 Calibration of electro-chemo-mechanical parameters for modeling solid state 


cathodes 


We present here a brief summary of the calibration of necessary parameters used to model 
the chemo-mechanical behavior of a LiCoO>-Li¡oGeP2512 composite cathode. In particu- 
lar, we performed calibrations to fit the Young’s Modulus, open-circuit voltage, and reac- 
tion constant for LCO. 

A sixth-order polynomial is used to capture the variation of Young’s modulus of Li, CoO» 
with Li content from E = 108.5 GPa for a pristine unlithiated state to E) = 252 GPa for a 


fully-lithiated state as reported in Wu and Zhang [140]. The result of our fitting is given by 


153 


E(x) = 136.932 — 84.93x* + 91.57% + 108.5 (A.3.1) 


where x denotes the stoichiometric amount of Li in Li CoO composition. 

Fig. A.3 (a) illustrates in red the analytical fit to the experimental Young’s modulus data 
reported in Wu and Zhang [140]. Also illustrated in Fig. A.3(b) is the variation in shear 
modulus with stoichiometric Li content employed in our simulations and compared against 
experimental data by Wu and Zhang [140]. Note that we did not fit for the variation in Shear 
modulus shown in Fig. A.3(b), rather it is produced from our fit for Young’s modulus and 


use of a constant Poisson’s ratio of y = 0.22. 
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Figure A.3: Analytical fit of (a) Young’s Modulus and (b) Shear Modulus for Li,CoO, 
active material with stochiometric Li amount - x to experimental data by Wu and Zhang 
[140]. Reproduced with permission from [17]. 


The activity coefficient of Li for an LiCoO» electrode is calibrated by fitting a seventh- 
order polynomial (c.f [115]) to the open-circuit potential data for Li,CoO, against a Li 
reference electrode reported in the work of Mizushima et al. [145]. Fig. A.4 illustrates in 
red the analytical fit to experimental data. Also illustrated is the finite-element simulated 
open-circuit potential curve (yellow curve), demonstrating consistent results with experi- 
mental data. 

The reaction constant, kọ is calibrated to charge-discharge curves for Li,CoO, by Zhang 
et al. [152]. We adjust ko such that energy dissipation during a full cycle (evaluated as 


the area inside the Voltage vs. SOC curve) matches with the area from experiments. A 
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Figure A.4: Comparison of analytical, seventh-order polynomial fit and the finite-element 
simulated open circuit potential curve against the experimental data by Mizushima et al. 
[145]. Reproduced with permission from [17]. 


representative figure of our calibration scheme against the experimental charge-discharge 
curves by Zhang et al. at a C-Rate of C/2.7 is illustrated in Fig. A.5. Here, analytical, finite- 
element simulated and experimental Voltage vs. SOC curves are included and illustrate the 


consistency of our calibration procedure. 
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Figure A.5: Calibration of reaction constant, ky against experimental data by Zhang et al. 
[152], demonstrating consistency across analytical, numerical and experimentally obtained 
Voltage vs SOC curves. Reproduced with permission from [17]. 


A.4 Modeling of interfacial damage in a LiCoO>-LGPS composite cathode for an 


all-Solid-State Battery. Additional simulations. 


We include here additional simulation results, which for conciseness, were not included in 


the main body of the thesis and are presented here for completeness. 
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Role of particle size and size distribution on electrochemical performance and mechanical 


degradation. 


In addition to the simulations shown in Sect. 3.3 investigating the role of particle size 
and size distribution at a volume fraction of dam = 30%, we performed an additional set 
of simulations at dam = 20%. The generated results are included here for completeness. 
Fig. A.6 shows the simulation RVEs, all of which have a constant active particle volume 
fraction of dam = 20%, while particle size distribution is allowed to vary. As in the 
main body of this manuscript, particle size and aspect ratio are seeded using a uniform 


distribution with prescribed lower and upper bounds. 
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Figure A.6: Simulated RVEs with a constant active particle volume fraction dam = 20%. 
The particle sizes are given a uniform distribution with lower and upper bounds in the range 
of: i) 2-6 um 11) 4-5 um and iii) 6-7 um. The aspect ratio lower and upper bounds are set 
to 0.9 and 1.5 respectively for all cases. Reproduced with permission from [17]. 


The results of these simulations are shown in Fig. A.7, where again we observe con- 
sistent results in electrochemical performance and damage evolution across the three sim- 
ulated RVEs with different particle size distributions and a constant active material volume 
fraction dam = 20%. The results are in agreement with those discussed in Sect. 3.3 in the 


main body of this thesis. 
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Figure A.7: Simulation results for composite cathode RVEs with constant active material 
volume fraction dam = 20% and varying particle size distribution. (a) Voltage vs. SOC 
behavior also compared to simulations with no damage. (b) Evolution of average Denem as 
a function of time for the first half-cycle. (c) Evolution of Voltage as a function of average 
Denem. Reproduced with permission from [17]. 
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APPENDIX B 
PRINCIPLE OF VIRTUAL WORK, MACRO- AND MICROFORCE BALANCE 
LAW 


We present here a detailed derivation of the macro- and microforce balance laws for the 
rate-like kinematical descriptors in our theory. We consider a list of generalized virtual 
velocity fields to be given by V = (9x,0F*, 0€, VÓ£, dd, Vód) and constrained through 
the kinematic relation (4.4.2) presented in Sect. 4.4.1. For any part P with outward unit 
normal n, of the reference body B, the internal and external power are given by 


OWent(P, V) = f 


tr(Mz) - OK dar + Jr. - OX du + | nó€ dag + J yôd dag 
ap P ap ap 


(A.1) 
5Win( PV) = ih (S°: ES + ESE + G - Vee + wd + C- Vid) de, 
P 


To derive the balance laws, we invoke the power balance requirement (c.f. Sect. 4.4.1) and 
demand that dWex(P, V) = dWin(P, V) for all generalized virtual velocities V. 

First, we derive the local macroforce balance. Let € = 0, VÓ£ = 0, ôd = 0, Vdd = 0 
such that the kinematic constraint (4.4.2) yields the relation 0F° = (Vdx)F~'F*. For this 


choice, the principle of virtual power gives, 


i: trng) - Íx dag + Jr. - ÔX dug = [ser dug = [ST Vox dug (A.2) 
ap P P P 


which by defining 


T, = S°F°* (A.3) 
and applying the divergence theorem on (A.2) leads to the macroscopic force balance 


Div Tk + bg = 0, and the traction boundary condition tz(nz) = Tne. (A.4) 
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As is standard, the Piola stress T, is related to the symmetric Cauchy stress T through 


Pis TE 


and for future use we may write S° = JTF*'. 


(A.5) 


Subsequently, we derive the local microforce balances. The first microforce balance 


associated with reaction is obtained by letting ôx = 0, 6d = 0, Vdd = 0, such that (4.4.2) 


yields, 


S°: JF* = S°: ( = MEEN!) de 
= ( -F9S*: (PEN) Je 
= ( — (JETTE): (hEN) )os 
= ( Es JAM: N’) s£ 


where we have defined the elastic Mandel stress as, 
WW def FET TE? 
For this choice, and accounting for (A.6), the virtual power balance then yields 


nócdaz = | (=J h(E)M":N'SE + ESE + G- VE) dup. 
OP 


P 


Applying the divergence theorem leads to 


| (—J°h(E)M®: NOE + ESE — DivG ôE) dup = 1 (7 — G- ng) 6&dag 


P oP 


(A.6) 


(A.7) 


(A.8) 


(A.9) 


Granted this relationship must hold for all P and all values of ô£, standard variational 
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arguments yield the second microforce balance, 


E-—J“h(£)M*:N' — DivG = 0 (A.10) 
along with the corresponding boundary condition 
(Mg) = G > ne. (A.11) 


Lastly, we derive the local microforce balance associated with damage. Let ôx = 0, 0€ = 0 


and Vó£ = 0. In a similar fashion, the virtual power balance statement then yields 


f yåd dag = [tose + €: Vód) dur (A.12) 
ap P 


Applying the divergence theorem on (A.12) gives, 


[cove — w)dd dug + f (y — ¢- ng) ôd dag = 0 (A.13) 


P ap 
which must hold for all P and all values of ôd, yielding the second microforce balance 
Div -w =0 (A.14) 
along with the corresponding boundary condition 
(Me) = Ç Mz (A.15) 


In summary, using the principle of virtual work we have derived one macroforce balance 
for the Piola stress, Tę as well as two corresponding microforce balances for the stresses 


{E,G,¢,a@} as outlined in Sect. 4.4.1. 
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